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IN THE UNITED STATES PATENT AND TRADEMARK OFFICE 
JHU/APL Docket No. 1443-NASA 

TITLE 

AUTONOMOUS SATELLITE NAVIGATION SYSTEM 

STATEMENT OF GOVERNMENTAL INTEREST 
The Government has rights in this invention pursuant to Contract No. NAS5-97271 
awarded by the National Aeronautics and Space Administration. 



BACKGROUND OF THE INVENTION 
FIELD OF THE INVENTION 

This invention relates to global positioning system based navigation techniques for 
developing position and velocity data required for autonomous operation of satellites 
15 incorporating on-board, event-based command systems. 

DISCUSSION OF THE RELATED ART 

Spacecraft tracking and navigation began when The Johns Hopkins University 
Applied Physics Laboratory determined Sputnik's orbital parameters as a function of the 

20 doppler shift of an on-board radio beacon. That was in 1957. Soon thereafter, that 

organization devised and developed the Navy Transit Satellite Navigation System which 
provided global navigation services for the military and civilian community for over 35 
years. The "TRANSIT" system, used a network of low earth orbit satellites and therefore it 
was unacceptable for most spacecraft navigation applications. 

25 A limited number of experimental systems developed by the Applied Physics 

Laboratory (APL), NASA and the U.S. Air Force provided semi-autonomous means for 
orbit determination and prediction. The experimental system provided ground-based orbit 
determination for spacecraft through the use of complex ground stations and coherent 
transponders that reside on the spacecraft to be tracked. The French National Space Agency 

30 has operated a similar system, DORIS, since the early 1 990s. However, the high cost of 
operating such systems has forced spacecraft operators to find other tracking methods. 
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Over the iasi three decades, APL and other organizations have demonstrated the 
feasibility of using GPS for positioning satellites and other high dynamic platforms. For 
over 20 years, the APL developed SATRACK system has utilized GPS translators and 
5 ground-based signal processing systems for trajectory reconstruction and guidance system 
evaluation of Navy Trident missiles. The first GPS-based navigation of a satellite occurred 
on Transat, an APL spacecraft launched in 1978 and operated for over 10 years. Transat was 
a Transit navigation satellite with an on-board GPS translator. 

Autonomous positioning of satellites using space bome GPS receivers was first 

10 demonstrated in the early 1980s when APL developed and flew four GPSPAC systems, and 
in the early to mid 1990s when NASA and the Jet Propulsion Laboratory used GPS receivers 
on the TOPEX/Poseidon spacecraft. More recently, other programs have adopted the use of 
GPS-based navigation systems for spacecraft. For instance, U.S. Patent No. 5,109,346 for 
"Autonomous Spacecraft Navigation System" issued to J. Wertz on April 28, 1992 describes 

15 a system using onboard observations of the earth, sun and moon to determine spacecraft 
attitude, instantaneous position and orbit based on multiple position estimates. Position and 
orbit data are derived by multiple deterministic solutions, including some that employ star 
sensors and gyros, and the multiple solutions are accumulated in a Kalman filter to provide 
continuous estimates of position and orbit for use when the sun or moon is not visible. 

20 Developing totally autonomous satellites results in large cost savings because of the 

elimination of ground support stations. In addition, ground supported satellites are subject 
to ongoing maintenance costs and are prone to serious malfunction if ground support 
stations are damaged or if control signal transmissions are interfered with. 

OBJECTIVES OF THE INVENTION 

25 A primary objective of the present invention is to provide a satellite for collecting 

data which utilizes on-board event-based command techniques relying on a navigation 
system which uses an extended Kalman filter to process data derived from a plurality of 
global positioning satellites. 

2 
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A further objective of the present invention is to provide an orbital platform for 
collecting data which utilizes an autonomous navigation system incorporating an application 
specific integrated circuit for processing signals from global positioning satellites in concert 
with an extended Kalman filter. 
5 Another objective is to provide an autonomous, on-board, event-based command 

system responsive to a GPS based navigation system for controlling data collection 
functions on-board low-earth-orbit and medium-earth-orbit satellites. 

A further objective of the invention is to provide an extended Kalman filter 
combined with an application specific integrated circuit for determining continuous real time 
10 earth-sun vector data, from a GPS constellation comprised of at least four satellites. 

A still further objective is to provide an autonomous navigation system whose 
primary input data is derived from at least four satellites of a GPS constellation. 

An objective of the invention is to provide an extended Kalman filter for interpreting 
data, from a GPS constellation comprised of at least four satellites and predicting ground 

15 contact event times therefrom. 

Another objective is to provide an application specific integrated circuit for 
processing data required for tracking global positioning system satellites. 

A further objective of the invention is to provide an extended Kalman filter 
combined with an application specific integrated circuit for interpreting data, from a GPS 
20 constellation comprised of at least four satellites, as continuous real-time orbital platform 
position and velocity. 

A still further objective of the invention is to provide an extended Kalman filter 
combined with an application specific integrated circuit for propagating state vectors to 
allow estimating future orbital platform position from a GPS constellation comprised of at 
25 least four satellites. 

SUMMARY OF THE INVENTION 
The primary purpose of the present invention is to provide an autonomous navigation 
system for satellites to be employed in NASA's Solar Terrestrial Probe Program for 
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measuring Thermosphcrc. Ionosphere, Mesosphcre, Energetics, and Dynamics (TIMED) 
The TIMED program is a remote atmospheric remote sensing mission to study the 
interaction of the Sun and Earth's atmosphere. This spacecraft serves as an orbital platform 
for four instruments that measure the basic state parameters and energy balance of the 

5 mesospheric, lower thermospheic, and ionospheric regions, i.e. MLTI regions of the 

atmosphere. The program investigates and documents the energetics of the MLTI region, 
i.e., its pressure, temperature, density, and wind structure and the relative importance of 
various sources and sinks of energy. To make the required measurements, a 625-km circular 
orbit inclined 74. 1 degrees with a nodal regression of 7200 per year is used. On-orbit 

10 instruments augmented by an array of ground-based instruments provide the basic 

measurements required by the mission. The instruments require precise knowledge of 
position and velocity of the orbital platform to perform their mission. This is accomplished 
by the autonomous navigation and time keeping system which provides position, velocity, 
time, and earth-sun vector data and notification of defined orbital events in real-time. It also 

15 generates predictions of events such as ground stations contacts. In addition, it generates 
orbital element sets which are down-linked for use by ground station antenna pointing 
systems. The navigation system contains two processors; a tracking processor which is 
responsible for controlling and interacting with GPS hardware in order to obtain raw 
tracking data and the navigation processor which is responsible for command and telemetry 

20 handling, determination of navigation solutions, generation of tracking acquisition aids and 
generation of the output data products. To simplify command control, the GPS based 
navigation system (GNS) is used to produce accurate estimates of position, velocity, and 
time to an event-based command system. 

The GNS uses SPS ranging signals broadcast from the constellation of GPS 

25 satellites. The SPS ranging signal, referred to as LI is Binary Phase Shift Key (BPSK) 
modulated. The modulation consists of two components that are modulo-2 summed: (1) a 
1.023 MHz Pseudo-Random Noise (PRN) code known as the coarse acquisition (CA) code 
and (2), a 50 Hz navigation message. The CA.code sequence repeats every 1 ms. The GNS 
receiver demodulates the received code from the LI carrier, and detects the time offset 

4 
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between the received and a locally generated replica of the code. The receiver also 
reconstructs the navigation message data. 

To compute TIMED spacecraft position, velocity, and time, the GNS determines the 
pseudorange to four or more GPS satellites/ 1 The^ro'pagation time to each GPS satellite is obtained 

5 by determining the difference between transmit and receipt times of the CA code. The pseudorange 
to each GPS satellite is computed by multiplying each propagation time measurement by the 
speed of light. 

The navigation message transmitted from each GPS satellite provides data which are 
required to support the position determination process. That includes information to 

10 determine satellite time of transmission, satellite position, satellite health, satellite clock 
correction, time transferred to UTC, and constellation status. 

Figure 1 illustrates the 660 kg TIMED spacecraft which is three-axis stabilized and 
NADIR pointing. The majority of the bus electronics are contained in an integrated 
electronics module (IEM). The IEM is a highly integrated system with major spacecraft bus 

Ys' "subsystems implemented on cards in a card cage to eliminate the more common 

self-contained modules or "black boxes". The subsystems that populate the IEM are on 
cards which have common mechanical interfaces based on the Stretch SEM-E form-factor to 
simplify the design process. Two IBM processors are employed, configured as single-string , 
redundant systems with a 1 553 bus connecting the two. 

20 A graphite epoxy optical bench is mounted on the spacecraft's zenith-facing surface to 

provide a thermally stable mounting surface for TIMED Doppler Interferometer (TIDI) 
telescopes and the attitude system's star trackers (required to satisfy TEDI attitude knowledge 
requirements). Two GPS Navigation System (GNS) antennas are mounted atop composite 
masts on the optical bench for performing on-orbit GPS attitude determination verification 

25 with a precision reference, or truth, data set. 

The TIMED mission and spacecraft top-level command and control requirements are 
satisfied by the GNS which is designed specifically for TIMED and implemented on 
Application-Specific Integrated Circuits (ASIC). Although designed for TIMED, the GNS 



5 
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will satisfy the functional requirements thai arc typically placed on spaceborne autonomous 
navigation systems for LEO or MEO hosi vehicles. 

The GNS includes redundant Standard Positioning Service (SPS) receiver systems 
with access to the GPS civilian ranging code called the coarse acquisition (CA) code that 
5 modulates the LI (1575.42 MHz) signal. The GNS is a state-of-the-art spaceborne system 
optimized for autonomous on-orbit operations. The GNS is radiation tolerant, has extensive 
command and telemetry capability, provides access to raw and intermediate data products, 
supports on-orbit software reprogramming, is designed to accommodate the large doppler signal, 
dynamic range resulting from orbital velocities, and implements robust signal acquisition, 
10 navigation, and orbit determination algorithms. 

BRIEF DESCRIPTION OF THE DRAWINGS 
Figure 1 is a simplified line drawing of the TIMED spacecraft illustrating the relative 
positioning of the major subassemblies. 
15 Figure 2 is a simplified block diagram of the TIMED GPS Navigation System 

(GNS). 

Figure 3 is a block diagram of the hardware implementation of the GNS. 
Figure 4 is a block diagram of one of the RF front ends. 

Figure 5 is a functional diagram of the baseband electronics subsystem identifying 
20 hardware and software components. 

Figure 6 is a simplified block diagram of the GPS tracking application specific 
integrated circuit (GTA). 

Figure 7 is a schematic diagram of one of the four input latch and I/Q generator 
channels of the GTA. 

25 Figure 8 is a functional block diagram of one of the 12 tracker channels of the GTA. 

Figure 9 is a block diagram of the carrier numerically controlled oscillator, NCO, 
which functions as a local oscillator for the final down conversion at the head end of a 
tracker channel. 



6 
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Figure 10 is a block diagram of ihc numerically controlled oscillator used lo clock 
the GPS coarse acquisition code. 

Figure 1 1 is a block diagram of the GPS coarse acquisition code generator. 
Figure 12 is a block diagram of the clock and interrupt generation unit. 
5 Figure 13 illustrates the time error conventions related to trajectory position. 

Figure 14 is a graphic depiction of Kalman crank interval subdivisions. 
Figure 15 is a graph of the MIC slew control. 

Figure 16 is a functional block diagram of the dual processor subsystem 
implementation of the GNS software. 
10 Figure 17 is a simplified hardware layout of the dual processor subsystem. 

DESCRIPTION OF THE PREFERRED EMBODIMENT 
The TIMED spacecraft major components include a system for sounding of the 
atmosphere using broadband emission radiometry to measure pressure, temperature and 
15 infrared cooling. This is the SABER instrumentation package. Located on the underside of 
SABER, Figure 1, is a Global Ultraviolet Imager (GUVI) for obtaining temperature profiles 
and auroral energy inputs. A pair of NADIR S-band antennas are located with the GUVI. 

A graphite epoxy optical bench is positioned on top of the SABER instrumentation 
package. The bench serves as a mounting base for the TIMED Doppler Interferometer 
20 (TIDI) profiler and four TIDI telescopes that constitute the TIMED Doppler Interferometer 
system which obtains wind and temperature profiles. A pair of star cameras, positioned on 
the bench, provide inputs to attitude star trackers which supply the attitude knowledge 
requirements of the TIDI. Also positioned on the optical bench is the Solar Extreme-UV 
Experiment (SEE) which measures solar x-ray, UV and FUV irradiance and a pair of S-band 
25 antennas. 

Power for the system is derived from a battery power supply augmented by two solar 

arrays. 

The primary electronics for the TIMED system are contained in an integrated 
electronics module (IEM). The system is redundant, being comprised of two IEMs, IEM #1 

7 



BNSDOCID: <WO_0063646A1_1_> 



WO 00/63646 



PCIYUS00/10197 



and IEM #2, see Figure 3. Each IBM includes a 2.5 CJbil SSR with command and telemetry 
interface on a command and data handling card. An RF telecommunications card includes 
an uplink receiver and down link transmitter. A DC/DC converter card, 94, contains a pair 
of voltage converters. Also contained in each IEM is a GPS navigation system, GNS. Each 
5 GNS includes a GPS receiver card, 10, and a processor card, 90, containing dual Mongoose 
V processors . 

The GPS navigation system provides: 1) estimated position and velocity of the 
satellite; 2) estimated UTC time and transfer to the command and data handling card; 3) 
estimated earth-sun vector data; 4) position-based event notifications; 5) position-based 
10 event predictions; 6) autonomous orbit determination; 7) data for non-GPS navigation; and 
8) support for on-orbit software up-loading. 

The GPS navigation system, i.e., GNS, operations and specifications are presented in 

Table 1 . 



15 Table 1. GNS Operations and Specifications. 

Parameter Operation or Specification 



Position (P) 

20 Velocity V) 

Coordinate frame 
Earth-Sun vector (S) 
Time of validity of PVS 
Time, output signal (T) 

25 Time, output data (T) 

Data rate of P VST 
Event notification (E) 



< 300 m (3 o), Geodetic latitude, longitude, height (<10 m 
predicted) 

< 25 cm/sec (3 o), east, north, up (<5 cm/sec predicted) 
ECEF-CTS (ECI-CIS for use by attitude system) 
0.06° (3 a) 

UTC 1-sec epochs 

100 ^sec(3 o)(< 1 ji sec predicted) 

CCSDS Unsegmented Time Code (CUC) Epoch: OH, Jan. 6, 

1980 

1 Hz 

Notify spacecraft of events, 5 sec uncertainty 



8 
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({veins: i erminator crossing: primary and backup ground coiuads: 

SAA; polar region 
Orbit determination Every 12 hours, predict forward 60 hours 



5 Data format: 



Event predictions 



NORAD 2-line element set and one state vector per ground 
contact 

Every 12 hours, predict forward 60 hours 



Non-GPS navigation 

10 

Autonomous operations 
Uploads accepted: 



Input voltage 
1 5 Power dissipation 

Physical configuration 
Radiation 



No accuracy reqs; errors grow w/atmospheric density 
uncertainty 

TTFF <30 mm, (cold start) 90% probable; <2 min (warm start) 
Solar flux; Geomag. Index; Polar wander; UTI-GPS time 
offset 

±15 Vdc and ±5 Vdc 

= 7W (1.5 W not including processors) 

Two 2-sided stretch SEM-E circuit boards + preamp + antenna 
5 krads (Si); latch-up immune (GTA & proc> 1 Mrad, RF 
IC: ;> 50 krads) 



20 Figure 2 is a simplified block diagram of the GNS. The receive antenna, 1 1, is 

located on the optical bench on the zenith-pointing surface of the orbital platform. The 
preamplifier, 12, consisting of pre-select filters and a low noise amplifier, is located just 
underneath the optical bench. The preamplifier generates the dominant component of 
system thermal noise power and conditions the received signals and noise for processing by 

25 the RF downconverter, 13. The downconverter, using a low-noise, 

temperature-compensated crystal oscillator (TCXO) as a frequency reference, downconverts 
L-band signals to the baseband and converts the resultant analog signal to a two-bit, 
three-level digital signal. 



9 
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Ai the heart of the electronics is a low-voltage. low-power, radiation-hardened, 
high-speed application specific integrated circuit (ASIC), 15. This GPS Tracking ASIC 
(GTA) implements all the GPS-specific digital hardware functions. 

The received GPS signals are converted to. baseband digital signals and input to the 
5 GTA from the downconverter, 13, via a baseband electronics interface board. The GTA, 
under software control, implements 12 independent GPS tracking channels, 16, provides 
timing and control functions, and incorporates a built-in-test capability. The device is 
designed to accept up to four downconverted inputs to support future attitude determination 
capabilities and can be daisy-chained to support over 72 tracking channels. It will support 

10 P/Y code operations with additional external hardware. 

The GNS's computer system, 90, comprises two Mongoose V radiation-hardened 
32-bit processors, the tracking processor (TP), 91, and the navigation processor (NP), 92, 
and associated memory and interface support circuitry. Each Mongoose V operates at 12 
MHz with no wait states for memory access. The two-processor approach is used to ensure 

15 adequate processor throughput. In an alternate embodiment, a single processor, such as a 
20-MHz Mongoose V processor is used. The tracking processor, 91, implements the GNS 
acquisition, tracking, and data preprocessing functions, and outputs its products to the 
navigation processor, 92. Optimized "lost in space" or "sky search" signal acquisition 
algorithms may be executed on the tracking processor. These algorithms are designed to 

20 minimize the cold-start time-to-first-fix and to minimize the likelihood of on-orbit GNS 

signal acquisition problems. The navigation processor, 92, implements numerous functions 
including the command, control, and telemetry processing, and the navigation and orbit 
determination algorithms including a 45-state extended Kalman filter, orbit propagator, 
event detector and predictor, and an orbit element set generator. 

25 The GNS is implemented redundantly on two identical TIMED Integrated 

Electronics Modules (IEM), each of which includes a GNS receiver card, 10, and a GNS 
dual-processor card, 90, see Figure 3. The receiver card, 10, and processor card, 90, 
interface to each other and to the rest of the IEM via a backplane bus, 93. Regulated DC 
power, 94, and discrete reset signals are input to each card and discrete analog telemetry' 

10 
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signals and a 1 PPS signal, steered to align with UTC, are output via the Cmd & Telemetry 
card, 95 . The interface to the IEM's command and data handling (C&.DH) system, 96, 
including input commands and all output data products, is provided through a 16-bit, 
2.5-MHz (40 Mbps) implementation of the PCI bus through the backplane connector. 
5 Six major elements comprise each GNS, they are: 1 ) a GNS micro-strip antenna, 1 1 ; 

2) an RF subsystem on a GNS receiver card, 10; 3) a baseband electronics subsystem on the 
GNS receiver card, 15; 4) a dual processor subsystem on an independent GNS dual 
processor card, 90; 5) an IEM backplane bus I/O network, 93; and 6) the system soft-ware, 
Figure 16. 

10 

GNS Antenna 

Each GNS antenna is a single-element microstrip patch bonded to a truncated, 
cone-shaped ground plane. The antenna is painted with a thermal control paint and coated 
with 1200 angstroms of Si0 2 to protect the paint from atomic oxygen erosion on orbit. The 

15 element is tuned, after compensation for the effects of the paint, to receive right-hand 

circular polarized (RHCP) signals at the GPS LI frequency of 1575.42 MHz with a nominal 
bandwidth of 5 MHz. In the preferred embodiment, the antenna is a modification of a ^ 
space-qualified design produced by Ball Aerospace and Technology Corp. The shaped 
ground plane provides excellent gain (nominally ^ 2.5-dBiC) at ^10° elevations. Antenna 

20 gain on the order of +5 dBiC is achieved at zenith. To eliminate obstruction of the GPS 

signals and to minimize multipath effects, the antenna is mounted on top of a graphite epoxy 
pedestal, which in turn, is mounted on the zenith optical bench, providing a full hemisphere 
of unobstructed coverage. 

25 RF Subsystem 

The RF subsystem, Figure 3, is comprised of a GNS Receiver Card, 10, for each 
GPS antenna, 11. Each receiver card, 10, includes four RF preamplifier modules, 12, which 
are space-qualified COTS devices with a preselect LI filter and a low-noise amplifier, 122 
of Figure 4. The filter is comprised of a 3-dB bandwidth at 30 MHz element, 1231 and a 

1 1 
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90-dB bandwidth at 250 MHz element, 123. It amplifies the band-passed raw received 
signal and applies it to an RF, L-band downconverter module, 13. 

The downconverter, 13, is based on a COTS triple-superheterodyne downconverting 
integrated circuit, such as a Plessey 6P2010, adapted for the GNS design. It is a plastic 
5 encapsulated microcircuit (PEM). There are eight downconverters, 13, on a GNS receiver 
card, 10, but only one is illustrated for simplicity as in the case of the eight RF preamplifiers, 
12. 

A PPL frequency synthesizer, 19, functions as the local oscillator for the 
downconverter. It is driven by the fundamental reference oscillator, 17, which is a 

10 high-stability, low-noise TCXO crystal oscillator. The oscillator's phase noise is less than 
45 dBc at 1 Hz offset from the carrier and -95 dBc at '10 Hz with stability of less than ± 1 
part per million (PPM) over all operating temperatures. The reference oscillator, 18, also 
generates the system clock for the GTA. A surface acoustic wave (SAW) filter, 18, in the 
second intermediate frequency (IF) circuit provides filtering of out-df-band signals and noise 

15 and establishes the system bandwidth. 

There are four RF receiver channels on each GNS receiver card, 10. Each receiver 
channel comprises two RF preamplifier modules, 12, and two RF downconverter modules, 
13, supplying a set of complementary I and Q signals. Thus the eight RF 
preamplifier/downconverter circuits provide sampled two-bit, three-level signals in four sets 

20 of complementary I-signals, compl and comp2, and complimentary Q-signals, compl and 
comp2, for digital processing. 

A block diagram of one of the RF front ends is presented by Figure 4 to illustrates 
how the GPS signal is downconverted and digitized in the preamplifiers, 12, and 
downconverters, 13. These modules produce a single pair of complimentary outputs as I- 

25 signals or Q-signals. After downconversion, the signal is buried in noise but it is extracted 
as the state of a three level quantizer, 14. The outputs of the eight quantizers, as I-signal or 
Q-signal compl and comp2, are applied to the baseband electronics module, 15 of Figure 3, 
as the four GPS receiver channel output signal pair sets, compl and comp2 for I and Q 
signals. The signals are applied to an associated one of the four input channels of the GPS 
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tracking application specific integrated circuit, 16, which is part of the baseband electronics 
system, 15. 

Baseband Electronics Subsystem 
5 The baseband electronics subsystem is comprised of the hardware and software 

elements illustrated in Figure 5. The software element is an extended Kalman filter process 
executed by microprocessors on the processor cards 90. The hardware is comprised of 
processor cards 90 and receiver cards 10 which are constructed around a GPS tracking 
application specific integrated circuit or GTA. 

10 Figure 6 is a simplified block diagram of the GTA. It is implemented in radiation 

hardened CMOS circuitry providing 12 independent tracker channels, 20, for acquiring and 
tracking the downconverted GPS signals from four receiver channels. The GPS 
downconverted signals are applied to the GTA input which comprises four identical IF 
input channels, each of which includes a data latch, 22, and an I/Q generation module, 23,, 

15 for adapting various input formats to the needs of the GTA circuitry and program routines. 
The I/Q generation modules supply quadrature and non-quadrature in-phase data to a data 
router, 21, which feeds the data to an available one of the 12 tracker channels, 20, as a 
function of channel availability and GPS satellite source. The GTA is relatively self 
contained in that it includes a time base, a clock, an interrupt generator, 24, and a 16-bit _ 

20 microprocessor memory-mapped interface, 25. Each of the 12 tracking channels, 20, 

incorporate full in-phase and quadrature (I&Q) tracking loops that correlate local replicas of 
the GPS PRN codes with the received signals. When acquiring and tracking GPS signals, 
numerically controlled oscillators (NCOs) are steered by the tracking processor's control 
loop filters to force alignment of the received PRN codes and the locally generated replicas 

25 clocked by the NCOs. Similarly, the reconstructed transmitted carrier signal is tracked by a 
carrier tracking loop. Each of the 12 channels also serve as a signal search engine when new 
GPS satellites come into view of the GNS antennas. 

The GTA handles a wide variety of sampled data inputs. The expected forms of 
sampled input to the GTA include both 1 bit and 1.5 bit versions of quadrature (I&Q) and 
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non-quadrature (i only) outputs from the downconverters. Quadrature inputs are desirable 
when the downconve'rsion process results in spectral folding due to a lower than nominal 
sampling rate. The losses of spectral folding are recovered when quadrature samples are 
used in a. single side-band (SSB) carrier mixing process which removes the image frequency 
5 in the signal band. If spectral folding does not occur in the downconversion process, the 
non-quadrature (I only) samples are adequate. 

Regardless of the GPS signal data type received from the RF downconverters, 13, it 
is applied to one of the four input channels to the GTA. Figure 7 is a schematic diagram of 
the input latch, 22, and I/Q generator, 23, circuits which form each of the four input 

10 channels. The latch circuit includes a pair of exclusive NOR gates, 3 1 and 32, which route 
the data in separate non-quadrature, I, and quadrature; Q, paths to enable the system to 
obtain an ideal set of quadrature sampled data. The I and Q input signals are each comprised 
of complimentary outputs compl and comp2. This approach is used because of its inherent 
immunity to phase imbalances which exists between carriers in the quadrature 

15 downconverters. Each exclusive NOR gate, 31 and 32, provides inputs to the I/Q 

multiplexers 33 and 34. The multiplexers, 33 and 34, are each responsive to both exclusive 
NOR gate outputs and the comp2 inputs in the I and Q channels. The compl I and Q 
channel inputs are applied to a third multiplexer, 35. The three multiplexers, 33, 34, and 35, 
form the channel's I/Q generation module, 23, which enables the channel to handle different 

20 input formats and convert them to an output which the remainder of the GTA can recognize. 
There are eight different modes of operation which the I/Q generator can handle. Each mode 
is listed and described in Table 2 in terms of the output produced for the reset of the GTA. 
This table includes the definition of the 1.5 bit output signal in terms of two actual bits. 
During the I-only operation (Modes 5-7), the zero bit of the Q sample is set. The I/Q 

25 generation module 1 .5 bit representations are presented in Table 3. 



14 



3NSDOCID: <WO 0063646A1 j_> 



WO 00/63646 



PCT/US00/10197 



Table 2. I/Q Generation Modes 
Mode Description 

0 1 & Q, 1 Bit, use compl input of both I & Q channels 

1 1 & Q, 1.5 bit, generate zeros from compl and comp2 inputs 
5 2 I & Q, 1.5 bit, comp2=zero bit function for both I & Q 

3 I & Q, 1 .5 bit, generate quadrature data from non-quadrature data 

4 I & Q, 1 bit, generate quadrature data from non-quadrature data 

5 I, 1 bit, use compl input from I channel only 

6 I, 1.5 bit, generate zero bit from compl and comp2 inputs 
10 7 I, 1.5 bit, comp2=zero bit function for I only 



Table 3. 1.5 Bit Representations 
Signal Value Sample (1) Sample (2) 

-10 0 
15 0 1 X (don't care) 

1 0 1 



The data router 21 acts as an intelligent distribution means for the outputs of the I/Q 
generators, 23. It provides an I sample and a Q sample input pair to one of the 12 tracker % 

20 channels, 20. As seen in Figure 8, the functional block diagram of one of the tracker 
channels, each tracker channel includes an I sample input, 41, and a Q sample input, 51 , 
from the data router. The I and Q sample channels are identical so only the I sample channel 
is discussed in detail. It is comprised of a suming network 43 which combines the outputs of 
two carrier mixing operations, 42 and 45. The first mixing operation, 42, multiplies the I 

25 sample input, 41, with a cosine function from the phase mapping generator 44. The other 
mixing operation, 45, multiplies the Q sample input 51 with a sine function from the phase 
mapping generator 44. The output of the summing network 43 is re-clocked with a flip-flop 
46 which provides parallel outputs to the three code mixers comprised of accumulators 47 
and associated latched registers 48 and an un-mixed output 49. 
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The Q sample channel is nearly identical to the I sample channel. It is identified in 
Figure 8 by reference designators 51, 52, 53, 55, 56, 57, 58 and 59 which provided functions 
similar to those provided- by. the circuit elements 4 1, 42, 43, 45, 46, 47, 48 and 49 of the I 
sample channel. The key difference is the summing operation 53 is indicating the difference 
between the mixing outputs 55 and 52. 

The outputs of the flip-flops 46 and 56 represent the final downconversion process 
which is achieved using an adjustable local oscillator frequency to accommodate the variable 
Doppler shift present in the GPS downlink. A numerically controlled oscillator, the carrier 
NCO, 54, generates a digitally responsive version of the local oscillator frequency required 
for the final downconversion process. A block diagram of the oscillator is presented in 
Figure 9. It uses a 24-bit full adder, 67, to provide the sequential logic which is clocked at 
the sample rate of the signal by latch 61 and kept in phase with sample by the phase and rate 
of phase change circuit 68. The three most significant bits, MSBs, of the adder provide the 
phase output, 62, to the mixing process in the sample channels via the phase mapping means 
44. The carrier NCO provides the system with an ability to vary local oscillator frequency. 
A register, 63, allows the system to read the current phase information. The system can set 
the NCO output to Os at anytime. The carrier NCO has a 32 bit counter, 64, with which it 
keeps track of the number of carrier cycles occurring during each measurement period and 
provides integrated Doppler estimates during GPS tracking. It's value is read by the 
microprocessor on each measurement interval. The carrier cycle counter and NCO phase are 
set to 0 by either an align_ time condition or an nco_car_clear condition. The controls are 
maintained by the microprocessor and actions occur only on the first positive edge of the 
MIC signal following a clear request. The 3 bit phase output of the carrier NCO, 54, is 
mapped to sine and cosine components according to Table 4. 
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Table 4. 


Carrier NCO Phase Mapping 






Carrier NCO Phase 


Cosine Value 


oine vaiue 




000 


+2 


+ 1 




001 


+ 1 


, ,/ 


5 


010 


-1 






on 


-2 


+1 




100 


-2 


-1 




101 


-1 


-2 




110 


+1 


-2 


10 


111 


+2 


-1 



Three fundamental signals act as inputs during the carrier mixing process in Figure 8. 
The three signals are the two quadrature sample data components, Isamp, 41 , and Qsamp, 
5 1, and three bits of phase from the carrier NCO, 54. The two quadrature carrier 

15 components, cosine and a sine, are generated by the phase mapper, 44, from the carrier NCO 
phase information according to the data presented in Table 4 . The carrier mixing process is 
similar to processing employed for single sideband modulation (SSB). This technique 
provides a large amount of image rejection. It is used to prevent image products from 
folding back down to the desired signal band and corrupting the signal. The fundamental 

20 mixing process is represented by the equations: 

/ = Isamp* cos+ Qsamp* sin 
O = Isamp* sin- Qsamp* cos 

The carrier mixing process is implemented with multiple-bit accuracy. The input 
samples, Isamp and Qsamp, have the values -1, 0, or +1. The sine and cosine components 
from the carrier NCO have the values -2, -1 , +1 , and +2. This results in mixing outputs 
25 which can take on the values 0, ± 1 , ± 2, ±3, and ±4. These values are represented by a 

four bit binar/ output defined in Table 5. The outputs are re-clocked at the sample clock rate 
for timing considerations. 
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Tabic 5. Mulii -Value Bil Definition 



5 



10 



Value 


13 i nary 


-4 


0000 


-3 


0001 


-2 


0010 


-1 


0011 


0 


1XXX 


1 


0100 


2 


0101 


3 


0110 


4 


0111 



15 The numerically controlled oscillator, NCO, which clocks the GPS course 

acquisition code generator circuit, the CA code NCO, 88, is similar to the carrier NCO, 54. 
A block diagram of the CA code NCO 88 is presented in Figure 10. To emphasize the 
similarity between the NCOs 54 and 88, common reference designators are used to indicate 
identical circuits in Figures 9 and 10. For instance, in Figure 10, the CA code NCO phase 

20 and rate of phase change is set by register, 68, which is identical to register 68 in Figure 9. It 
is also controlled by the microprocessor bus as in the carrier NCO. The phase increment 
register, chx_code_phincr, is used to set the phase offset incremented to the current phase at 
every sample clock time as in the carrier NCO. The current phase clock of the CA code 
NCO, 88, is read at the current phase register, 63, by performing a microprocessor read from 

25 the address chxjcodenco jphase as with the current phase register 63 of the carrier NCO, 54,. 
Verification of correct CA code NCO operation is determined by directing the 
microprocessor to clear the current phase prior to the next phase addition by setting the 
appropriate bit in the register nco_code_clear. This clear operation occurs at the next 
measurement interval clock, MIC, cycle. 

IS 
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The GPS coarse acquisition code numerically controlled oscillator, CA code NCO, 
88 of Figure 8, provides an output to the GPS course acquisition code generator 70. A block 
diagram of generator 70 is illustrated in Figure 1 1 . It incorporates two shift registers, 71 and 
72, which generate the PRN codes that are summed together by the modulo-2 summing 

5 means 73 to produce the desired GPS coarse acquisition code. Each shift register is clocked 
on the negative edge of the CA code NCO fchip to ensure that the code generator is 
initialized at a specific time so the current phase of the code may be derived from a code 
phase counter. The code phase counter, 74, is a 10-bit counter that resets to a value of zero 
at every epoch and increments by one on every chip. Critical timing is achieved with the use 

10 of the align Jime signal form the microprocessor along with the MIC signal. Different codes 
corresponding to each GPS space vehicle are generated by verying the initial phase of the 
lower shift register. It is possible to start generation of a GPS course acquisition code at a 
point besides its epoch by initializing the top register to something other than all Is. 

The next step in processing consists of mixing the signals, i.e. the outputs of the 

15 differentiators 46 and 56, with a locally generated replica of the code used during the data 
modulation produced by CA code generator 70 . This process occurs simultaneously in six 
differential mixing circuits, 81 through 86, for every tracker. The purpose is to lock 
frequency and phase of the'eode present on the received signal with a locally generated code. 
Both the I and Q outputs of the carrier mixing are multiplied by three time-shifted versions \ 

20 of the local code produced by the 3-bit shift register 87. This results in a total of six code 
mixing and correlation operations per tracker. Early, prompt, and late time-shifted versions 
of the code are used because their 1/2 chip delay offsets provide good code acquisition and 
tracking performance. 

The code mixing process is illustrated in the functional block diagram of a single 

25 tracker in Figure 8. The components of the I and Q carrier mixing products are clocked at 
the input sample rate and the code to be mixed with them is clocked at twice the CA code 
numerically controlled oscillator rate, i.e., CA code NCO, 88, rate . Because the CA code 
NCO is clocked at the sample rate, data transitions for all of the involve signals occur on the 
same boundaries. Both the 1 and Q carrier mixing products are adjusted so they are 

19 



BNSDOCID: <WO 0063646A1J_> 



WO 00/63646 



PCT/USOO/10197 



represented by the values 0, ± 1, ± 2, ±3, and ±4. Using the code values ± 1, the output of 
the code mixers, 81 through 86, take on the same range of values as the outputs of the carrier 
mixing process. 

After the code mixing process, the GPS signal outputs IE, IP, IL, QE, QP, and QL 
5 are loaded in the registers forming the appropriate accumulators, 47 and 57. Because the 
accumulator can be implemented easier as an up-counter instead of an up/down-counter, the 
representation of the mixer outputs is changed by adding an offset of +4 to the mixer 
outputs. The new values represent the increment value corresponding to each possible code 
mixing output. The process is repeated for subsequent values of the products coming out of 
10 each code mixer. This represents the correlation operation of a GPS receiver. See Table 6. 

Table 6. Relationship of Code Mixing Output to Counter Increment 





Code Mixing Value 


Code Mixing Bindery 


Counter Increment Value 


Counter Increment Binary 




-4 


0000 
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0000 


15 


-3 


0001 
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0001 




-2 


0010 
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0010 




-1 


0011 


3 


001 1 




0 


1XXX 


4 


0100 




1 


0100 


5 


0101 


20 


. 2 


0101 


6 


0110 




3 


0110 


7 


0111 




4 


0111 


8 


1000 



The six accumulator values from the code mixing outputs of each tracker are 
normally latched into a register 48 or 58, at every internally generated code epoch. 
25 However, to reduce the load on the microprocessor, each tracker has the capability of 
latching the data at a multiple number of epochs. Each tracker can be program so that 
anywhere from 1 to 20 epochs may occur between latching of the accumulator contents. 

The timebase, clock, and interrupt generator, 24, create all timing and interrupt 
signals for the GTA. Certain portions of the clock generator work specifically with the 
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TIMED frequency plan while other pans of the generator are adapted to provide different 
frequency and timing requirements. Figure 12 is a block diagram of the timing and interrupt 
generator. The generator is normally driven by the same sample clock that the rest of the 
GPA uses for data processing. 

5 A common clock or timing signal is required to efficiently allow the microprocessor 

to download the latched accumulator values between epoch states. This download signal, 
the accumulator or interval clock (AIC), occurs at a 1.499812 kHz rate or a period of 666.75 
microseconds. It is used to generate an interrupt for the microprocessor. 

The measurement interval clock (MIC) is a timing signal used to notify the 

10 microprocessor when measurement readings are to be be made. The MIC is used as a 
reference time for aligning other signals such as the carrier and CA code NCOs and code 
generator outputs. The MIC has four different, selectable rates, they are 100 Hz, 10 Hz, 1 
Hz and 0.1 Hz. During applications requiring critical timing measurements, a 1 PPS signal 
is used as the MIC. When using one of the internal MIC rates, this signal is a non-square 

15 waveform with a positive pulse width of two AIC clock periods, or 133334 microseconds. 

The time base generation circuit provides a 1 pulse-per-second (pps) UTC clock 
reference that is used to control the TIMED spacecraft timing system. Under control of a 
tracking loop in the navigation processor's extended Kalman filter, a numerically controlled 
oscillator is continually controlled such that it outputs a 1-PPS signal that is steered to align 

20 it with UTCs 1 -second epochs. This is accomplished by using a 25-bit synchronous divider 
(counter) which is clocked by the sample jzlk signal. When the counter reaches a value set 
by the microprocessor, the 1 PPS signal is set high and the counter resets to zero. The 1 PPS 
signal is then set low halfway through the next count to create a clock with a 50% duty 
cycle. This clock is capable of being aligned with UTC time by modifying the 25 bit 

25 maximum count value. Initially the 1 PPS divider is aligned with the MIC signal. This is 
achieved whenever an alignjtime signal occurs. 

The time offset between the MIC and a 1 PPS signals must be measurable. The 
measurement is obtained with a second counter that determines the time interval between the 
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1 PPS rising edge and the next MIC rising edge. Two registers hold the number of 
samplcjclk cycles which occurred between the two signals. 

The GTA is implemented on a radiation-hardened HR2300 series CMOS gate array 
produced by Honeywell Solid State Electronics Center. The GTA uses over 200,000 gates, 
5 operates on 3.3 V for internal circuits and 5 V for input/output, and dissipates = 200 mW.. 

Receiver Timing 

Receiver timing is unique for TIMED. The receiver contains a fundamental 
oscillator which drives a group of frequency dividers, one of which is coupled into the 
10 receiver clock to produce a function of true time which is uncontrolled or open loop time. 
Controlled time is created by a control process that drives the system time bias to zero. The 
process is modeled as a continuous time adjustment to frequency error which converts an 
uncontrolled time function to a controlled time function of true time for use by the Kalman 
filter. 

15 Correct modeling of the TIMED orbital dynamics and the GPS measurements is key 

to achieving good navigation performance. This includes having adequate gravity and drag 
models and correct modeling of the measurements derived from the GPS receiver carrier and 
code tracking loops. Modeling the effects of the receiver clock error on the measurements is 
critical. A unique feature of the GNS receiver is the control of the sampling times of the 

20 tracking loops to align with Universal Coordinated Time (UTC) time tics. 

Figure 4 illustrates how the GPS signal is downconverted and digitized by 
preamplifier 12 and downconverter 13. These circuit uses a Plessy RP front end and the 
GPS\ LI signal is downconverted from 1575.42 Mhz to 4.3 MHz. This technology leaves 
the signal buried in noise. The final output of the module is the state of a three level 

25 quantizer, 14. 

The outputs of the quantizer are processed by an application-specific integrated 
circuit, as shown in Figure 6. As previously stated, the ASIC chip includes twelve 
replications of a GPS tracker, so up to twelve GPS satellites can be tracked. The signals 
are sampled and distributed to the 12 trackers by four receiver channels. In the trackers, this 
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signal samples are digitally mixed wiih a local carrier and code. The sampling accomplishes 
an additional downconversion of the LI signal to -1 .4054 Mhz. The negative sign arises 
because the signal is oversampled and downconverted through zero frequency. This results 
in a reversal of sign in the Doppler shift. After lockup, the local carrier and code are in 

5 phase with the received carrier and code. The GPS message data is retrieved inside the 
carrier tracking loop as the pseudo-range and integrated carrier. The data is manipulated by 
an extended Kalman filter implemented by software in the mongoose microprocessors. 

The software used by the system is presented by way of mathematical functions to 
avoid hardware and programing language ambiguities. To simplify the presentation, the 

10 analysis considers the pseudo-range and integrated carrier produced by the code epoch 
accumulator and the carrier phase accumulator. using the output of a single exemplary 
tracker channel. These accumulators are sampled at interrupts generated by the 
Measurement Interrupt Clock (MIC). The timing of this interrupt is controlled by the 
Kalman filter so that at steady state, it occurs as close to the UTC time tics as possible. 

15 Since this is different than what normally occurs in a Kalman processed GPS evaluation, a 
thorough discussion of the receiver timing is presented. Receiver timing is developed by the 
GNS software which implements the following equations, where: 

Vectors/matrices are represented as lower/upper case {^/j^ • The transpose of a vector or 

matrix is denoted as x* . For quantities that depend on a particular tracking channel, the 
20 tracking channel superscript v is used. Thus a v is the pseudo-range associated with 
tracking channel v . Quantities associated with GPS satellites are identified by their 
tracking channel index (as opposed to their PRN or SV numbers). Thus x v (r) is the 
position implied by the GPS message ephemeris for the satellite being tracked by channel 
v . Estimates are denoted as x . Reference receiver time is denoted by / while true time is 
25 denoted by x . The reference times at the MIC interrupt event are denoted by t k . 

Sometimes t k is used to denote a general Kalman filter update time. (Note that these form a 
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subset of the MIC times since a measurement is processed every m' h MIC interrupt.) The 
apriori/apostcriori estimates of quantities at t k are denoted as x k/k _ , /x k/k . Unless 

otherwise stated, all vector quantities are expressed in the Conventional Inertial System 
(CIS). 

5 Basically, the system uses an imprecise clock. To develop error equations for a 

navigation system which contains an imprecise clock, a precise definition of error is 
required. This is accomplished in the software by defining the asynchronous errors as: 

S a x(t)^x(t)-x re/ (t) 

where x(-) is any function of time, x{t) is its true value, and x re/ (t) is its reference value 

10 according to the GNS system. Similarly, the synchronous errors are defined as: 

6 s x{t)=Mt))-x r «{t) 

where x(t) is a function that takes the nominal time t of the GNS system into the 
corresponding true time. This situation is illustrated graphically in Figure 13, where x{t)>t 
is assumed. If the difference x(t)-t is small, and the difference between the true trajectory 
1 5 and the reference trajectory is small, the following relations hold between the two error 
conventions: 

sAt)"SAt)*&Mt)-t) 0 ) 

where x(t) = C %j t 

The synchronous error is the error in x re/ (tj ignoring the time tagging errors of the GNS 
20 clock, while the asynchronous error takes these errors into account. If the difference 

between the indicated and true time to be the state is defined as T=t-x(t), the appropriate 
transformation between asynchronous and synchronous errors is: 
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S s x(r) 
T{t) 



1 - xit) 
0 1 



Tit) J 



and the associated suitably partition covariance is: 



cov 



(5 s x) cov(j 5 x 5 r) 
[symm] cov(T) 



I - x(t) 
0 1 



cov(^x) cov^x, t) 
(symm) cov(T) 



I - x{t) 
0 1 



For most sections in the following development, the synchronous error convention will serve 
5 as the default. It will be written without a subscript. 

For the uncontrolled time receiver reference times, Let r represent the independent 
variable true time. The receiver contains a fundamental oscillator with instantaneous 

frequency / 0 (r) ! . (For the TIMED GNS receiver using a Plessy RP Front End, f Q is 

nominally 10 MHz.) It drives a group of frequency dividers, one of which is coupled into 
10 the receiver clock. The output of this frequency divider as a function of true time is: 



sin 



K L 



r cw 



J/ 



= sin 



(f(r)/*,) 



where k c is the (rational) divide-down ratio. (For TIMED the system clock is coupled 
directly into the fundamental oscillator, so that k c = 1 ) The true time at which the clock- 
was set is r se( , the phase (in radians) of the fundamental oscillator is 
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(p{r) = 2/7 J" f Q (c)dc + <f> sel , and 0 sel e [0,2;r] is an arbitrary phase of the 

fundamental oscillator at r el . The number of complete cycles executed by this sinusoid 
since T set is N^J^t) = /loor[(<p(r) - 0 sel )/2 xk c ]. If / 0 rc/ cycles per second are 
being executed by the fundamental oscillator, the estimate of time is: 

v KM cycled ) 

* («■)= J77f + '«, 

/o 



= A e -y7oor[(^(r)-^ / )/2^ e ] ^ 

0 



where t sei is the value placed in the clock at r $et . Dispensing with the floor function since 
all times of interest (the sampling times) occur when the oscillator has executed an integral 
number of cycles (i.e. when {q>{r) - 0 sef )/2x is an integral multiple of k c ). 
Then: 



rref 
JO 

r /J j) 
0 



10 



(2) 



+ t 



set 
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where the fractional frequency error is defined by: £\T ) = 7^7 

/o " 

This is the uncontrolled or open loop time. If f wl = r wl and / 0 (^) = /J," 7 there will be no 

error and (r) = T . Since this is not the case, there will be an error defined by: 

r(r)=/*(r)-r 

= t sel - r se! - )e{t)dt W 

'sci 

5 Note that t* (r) is the function that takes true time into the simultaneous uncontrolled 
reference time. Note also that the divide-down ratio k c does not appear. This function is 
inverted to obtain the true time, r(t*) , as a function of uncontrolled reference time. This 
function is defined implicitly by: 

<>-) 

= = (*■('•)- O- Hfat- (5) 

10 The instantaneous phase of the fundamental oscillator is recovered from t*(r) by solving 
equation (2) for <p{t) to get: 

= 2xf 0 r «([r= T (r )]-/„, = j sel ) 
= 2*f«{T = r(r))+^ 
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where <j) sct has been modified to include t set , because only <p(r) is significant, modulo 2n 

and {t sct Inf™* ) mod(2/r) can be absorbed in 0 sc/ . 

The uncontrolled time thus far defined is used to derive controlled time by driving 
the system time bias to zero. A quantized and control in a discrete time system is modeled 
5 as a continuous time adjustment to the frequency error. Applying the control process to 
equation (2) produces: 

t(r)=t set + (r-r sel )- )[s(f)-4f)ty (6) 



where /(r) is the system reference time as a function of true time, including the control 

10 effects, and u(r) is the continuous time process used to represent the control. Its 

dimensions are seconds/second (the same as fractional frequency error). Note that t(r) is 
the function that takes true time into simultaneous controlled reference time. This function 
is inverted to obtain the true time, r{t) , as a function of controlled reference time. The 
inverse function is defined implicitly by: 



15 t = t m + 



(rM-rJ- J[*(fMf)]rff P) 



Defining the (controlled) clock bias to be that quantity which when subtracted from 
the (controlled) reference time gives the corresponding true time: 

7tr)=/(r)-r.. 

(8) 



= t sel 



- fl*G?M<f)]4f 
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The rate of change of T{t) with respect to true time is then: 

dT 



dr 



- -s(r)+u(r) (9) 



(10) 



The Kalman filter runs according to receiver time, therefor the derivative of 
T(t) = T(r(t)) with respect to t . By the chain rule, this is 

dT{t) _ dTjr(t)) dr{t) 
5 dt ~ dr dt 

f(t) = dr/dt is derived by differentiating equation (7) with respect to t using the 
Leibniz rule which produces: 

or 

10 r{t)=~ r-T p-r (11) 

where £-(r(z")) and = w(r(r)) . Substituting equations (1 1) and (9) into 

equation (10) yields: 

dT{t) _ ijt) - e{t) 
dt \+u(t)-£(t) (12) 

The variables r and / represent the true time and controlled reference time independent 
15 variables as well and the functions r{t) and its inverse t(r) . The distinction is clear from 
the context. When appropriate, superscripts identify the system to which the time functions 
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pertain. Thus, a reference time from the Guidance and Navigation System (GNS) leads to 
the corresponding true time: 

r(/)=/-7T(/) (13) 
By a notational slight-of-hand, as in equation (13 ), functions may be accessed 
5 according to receiver time instead of true time. This is accomplished by writing T(t) to 

mean T(r(t)) . The GNS time / represents the number of seconds in the GPS standard 

epoch. It is stored in two variables so that precision will not be lost: 

t k = MIC_time k + t m 

where MLCjime k is the number of seconds since GNS turn-on and t mx is an integer that 

10 brings the sum MICjime k + t inl close to the desired value when time is initially set. 

MIC slew control, u{t) , is a continuous time representation of the digital process 

used to control the timing of the MIC interrupt in the GNS receiver so that it aligns with 
UTC tic boundaries. The GNS receiver continually estimates the difference between its 
internal clock and GPS time. The GPS time differs from UTC by an integer; therefore the 
15 GPS time tics are coincident with UTC tics. In the actual digital implementation, the MIC 
interval counter is strobed by f Q /k c .An interrupt is generated and the counter resets 

when the counter reaches a value written into a memory location on the GPS Tracking ASIC 
(GTA). The desired number of clock cycles between MIC interrupts is known as the MIC 
register count and is determined from Kalman filter estimates. 

20 The objectives of the digital control are to compensate for the time bias error and the 

fractional frequency errors, and to bring the MIC interrupt into coincidence with the UTC 
tics. The control will be applied across the whole Kalman update interval consisting of/?? 
MIC interrupts so the difference between the true and nominal MIC interval times can be 
minimized. As described earlier, the nominal time since the GPS standard epoch is 

25 constructed as: 

h = MICjime k + ( in{ 
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The corresponding estimate of true time after a measurement update is found by removing 
the bias from the reference time: 

*kfk = MICjime k = / int - f k/k (14) 

It takes far less time to eliminate the clock bias and to align the MIC with UTC tic's 
5 if the clock bias is removed by absorbing it in t m : 

'fat «= rOUnd Vkl k) 

i- \ (,5) 

T klk <= T k/ k - round \T k/ k ] 

where round {x} is the closest integer to x , i.e. 

round {x} = arg min|z-x| (16) 

where Z are the integers. It is clear that this redefinition leaves the estimate of true time 
10 given by equation (14) unchanged. This redefinition is done after every measurement 

update. If the MIC slew is operating in steady state, the time bias will be maintained at a ' 
very small value and this redefinition will not do anything. The absorption of the integral 
part of the clock bias will only happen after the first measurement update or if the MIC slew 
is deactivated. 

15 Assuming x k denotes k th sampling time, then the k ih MIC interrupt is triggered 

when the GTA MIC counter (which is strobed at the frequency f 0 /k c ) counts up to a 
prescribed value, n{k) . The GTA MIC counter then resets. The true MIC interrupt times 
(r k ,k = 1,2,3. . . ) are thus defined by: 
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where r 0 is the true turn on time. 

The Kalman update interval contains m MIC intervals and a Kalman crank interval 

is: 

1 /»r k+m 

j- \;;-fii)d t - 2» 

5 The update interval is divided into subsections of length p and q = m - /? , as in 

Figure 14. 

For the first p MIC intervals, the measurement update calculations at t k are still 
underway and a default MIC interval count n d is used. The value for the default MIC 

10 interval count is based on the previously estimated fractional frequency error so that the true 
MIC interval times are as close to the nominal as possible: (If this is the first Kalman crank, 
the estimate comes from preflight calibration of the oscillator.) 

n d = ™w{a m/c ^ 

For the remainder of the Kalman crank interval, the MIC interval counts are calculated to 
15 minimize the time bias estimate predicted at the end of the crank interval: 

~r( T k +m " r *)= P" d + X (17) 

assuming f 0 (g) is slowly varying. Using the above definitions for the receiver time bias 
and frequency errors: 

= h ~ T{t k ) 

20 then 
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T k+m ~ t k+rn rfj k tm) 



and 



/ 0 = /o rC/ (l-^J) 
Substituting these values into equation (17) yields: 

f*{\-e(t k )) 



\ k+m 

~ T(t k+m )- [t k - T(t k )])= pn d + £ «(/) 



/"= k + p+ 1 



Let A M!C = t k+} - ^ be the nominal MIC interval in seconds (= lsec) . Then 



Minimizing to find the integer sequence n\i) which makes the clock bias at the next 
10 crank as small as possible: 



k + m 



Note that this optimization problem defines only the sum of the sequence S = ^ n(i) 



i=k+ p+] 



and can be expressed in the form: 
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S' = arg minj^ - SB\ 

= arg min\B(A/B- s)\ (18) 
= arg mm\A/B-S\ 



where: 



A= m-A wc + T(t k )- pn d 



l/o re/ (l-^))J 



B = 



/."(MO) 



5 By definition of the round function, the minimum of equation (18) is attained if : 
S' = round(A/B) 

V(M0). 



= round 



= round 



//^f 1 - / 4/ A J / X 

7 \m • A M/c + r t/ A j 



(19) 
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The last line is used for computation because the causal estimates are substituted for 
the true quantities z{j k ) and T{j k ) s which are unknown. 

The optimal sum is generate by the following sequence which is preferred over other 
sequences because it produces a minimum variance from the optimal mean count value 

5 S* jq . The last q MIC intervals in the Kaiman crank are divided into two groups, one 
with q f intervals using n{i) = floor^S* jq) as the MIC register count and one with 



q c = q - q f intervals with n{i) = floor(S* /q)+ 1 as the count: (20) 

n f = floor{s*/q) /e|/:+/?+l, k+p+2, ...k+p+q f ^ 
n c = floor(s* /q)+\ i e + p+ q f + 1, k+p+q f + 2, ...k+ wj 



»(/) = 



10 The sum of the sequence is: 

k+m 

E«W= q f n f + q c n c 

= q / -n / + (q-q / )(n / + \) 



i=k+p+) 



The integer q j that generates the minimum variance sequence that satisfies 

I< 

/ = fc + 1 



^ n{i) = 5* is: *S" = <3y • ^ + ( ^ - + 1 j which yields: 
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q f = q(n f + \)-S m 
which is greater than zero, provided q ,S* > 1 . Then q c is 

= S* -q n f 

During the period between the first two Kalman cranks, the MIC slew has to absorb 
5 time errors up to 1/2 second in magnitude. The MIC is disseminated to the spacecraft as a 1- 
PPS timing source. ThereTor, .the time between MIC interrupts cannot be shorter than 

0.9999 seconds. ( A MIC = 1 sec. ) To ensure that this constraint is met in the presence of 

a 100 nanosecond time quantization on the GTA and clock estimation uncertainties, the 
GNS uses a tighter constraint of 0.9999 + lOOx 1CT 9 + 3a T seconds, where a T is the 

10 estimation uncertainty of the receiver clock in seconds. Since A wc is expected to be 1 
second and a T is expected to be around 100 nanoseconds, S = 9.96 x 10~ 5 isused. If the 
shortest true time interval between MIC interrupts is less than A A/yc 0 - S) , the system 

redefines the MIC counts for MIC intervals k + p + 2 to k + m so that the true time 
interval between MIC interrupts is as close to A mc as possible and lengthens the 
15 k + p + \ th MIC interval so that the total length of the Kalman crank is increased by an 
integer number of seconds so that the length of the crank is raA wc + s seconds. 

The MIC interrupt has been brought into coincidence with UTC boundaries at the 
true time of the next Kalman crank, this is still the case if an integer number of true seconds 
are added to the Kalman crank interval. The rate of change of the time bias estimate in the 
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Kalman filter is computed based on the actual time sequence of MIC counts, assuring that 
the estimate of the time bias will evolve correctly. 

The best estimate of the shortest true time interval resulting from the calculations this 

far is n f k c j{^f^ cf {\ - £ k/k . The minimum MIC interval constraint will be violated 
5 (or\q f intervals) if: 



/o rc/ (i-^) 



< A„ /C (l- S) 



(23) 



If the foregoing is true, the system redefines the true MIC interval for MIC interrupts 
k + p + 2 to k + m to be as close to the nominal MIC interval as possible: 



k+ p+ 2, k + m 



(24) 



10 where 



/7 </.«/ = 



The foregoing is the value for the default MIC count that is used for the next Kalman crank 
interval. 

Next, the counts for the k + p\ Y h interval are defined so that the whole Kalman 

15 crank interval is lengthened by s seconds, according to the estimates of the clock parameters. 
A lengthening of s seconds means that instead of steering the clock bias to zero, the count is 
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selected for this new interval so that the clock bias 7l(/ jt+OT ) is as close to - s as possible. 
/?(/: + p + l) is selected from the integers to minimize: 



+ 5 = 



+ X "0") 



+ s 



k+m 



A + s- B X n U) 



i=k+p+\ 



A+s-B[n{k+p+ l)+{q-l)n dtnxt ]\ 
A'- Bn(k+ p+l)\ 



where A and B are as defined in equation (18) and A' - A + s- Byq - \)n d . The 
5 optimal integer value for n[k + p+ l) is then: 
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n(k+ p+\)= round { A'/ b) 

A + s- B(q- \)n d ^ 



= round^ 
= round \ 

= round 
- round 



B 



A+s 
B 



fo e/ ^-s k/k )(m-A MJC+ f k/k +s) 



~ P" d -{9- lK,« 



The integer 5 is chosen so that the true-time length of MIC interval is greater 

than A MIC {l- S) . This is accomplished by setting A r | > A MIC {} - S) 

or: 

n{k + p + \)>A wc (\-S) f °( ] ~^ (25) 
Since round(x) > x - 1/2 for x > 0 : + p+ l) > "'(5) where: 
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n'(s) = 



fo e/ ( 1 ~ e k , k ){m ■ A Aflc + f k/k +s) 

_ _ p „ d _ ^ _ _ 1/ 



k. 



If the magnitude of the oscillator fractional frequency error is bounded by \£ k | < /? , 
s is selected so that 

n(k +P+ l)> „is) > A wc (l - ^) /o " /( ,' + /?) > A mc (\ - S) ^ZlA 



5 or 



S > 



/o re/ (l-4/*) 



MIC 



+ pn d 



so 



^ = ceil 



/o re/ (i-4/*) 



7o re/ (i^) 



(<? - 0^.nxl + P n d ~ 171 • A MIC ~ 



k/k 



(26) 



The upper bound for 5 is defined by: 
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s < ceih 



'= ceib 



< ceib 



MIC 



f£!SliA {(l . sMq . l)+p)+ i. 



-m-A MIC - T k/k 



MIC 



(m- S) + 



2f«{\-5)_ 



-™ A MIC- T k/k 



m 



2^A 



MIC 



2f 0 re/ ) 



Now 



l k/k 



< sines the whole value part of the clock bias is absorbed in t mx (see 



equation 15), so s < 1 if: 



m ( n k c ) 1 



5 This implies 



n ef - mk c 



10 



With default values of m = 1 80, A M1C = 1 and f 0 re/ = 1 0 6 ; fi < 1.3867 x 10~ 3 

guarantees 5 < 1 . This is a very loose requirement for the oscillator fractional frequency 
error (and its estimate) since the expected flight oscillator will be bounded by 30ppm. 

The the following equations are used for calculation of the MIC interval counts for a 
Kalman crank interval. The default control used for the first p interval is: 



(27) ■ 
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then 



( 



B = 



= ™&Af/C +T k/k ~ P n d 



S' 



n 



f 



= round(A/B) 
= floor(S' /q) 



(28) 



and 



n c = rij + 1 



(29) 



If the minimum MIC constraint is satisfied (as it is in steady- state operation), i.e.: 

>k MIC {\-8) (30) 



n f k c 



using the baseline rules: 



«(/) = 



n d i e \k + 1, k + /?j 

n f i e{k+ p + 1, k + p+ #y j 

/7 c / g ^k + p+ q f + 1, £+/?-fmj 



(31) 
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where q f = q^n f + lj - 5" . If the minimum MIC interval constraint is violated, the 



following is calculated: 



n d** = found 



s = ceil 



/o re/ (l- V) 
A' =A+s-B{g-l)^ 



MIC 



(q - + P n d - m ■ A Mic - 



f (32) 



5 using the rule: 



n{i) = 



n b i e \k + 1, k + 

row«^{^7^} ze {/: + /?+ l} 

n nnxt i e{k + p+2, k + p+ w} 



(33) 



The average rate of change of the time bias over the Kalman update interval is: 



m> A 



A//C 
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Which is manipulated to: 



/; g/ (i- y) 

K 




to obtain: 




(34) 



•A w/c -/o re/ (l-^ A ) 



5 



This is the (constant) deterministic control value that is applied across the receiver 



time interval t t 



k+m 



when propagating the estimate of receiver time bias in equation 



(12). 



To illustrate the operation of the MIC slew control, consider first a case in which the 
MIC interval is 1 second, the Kalman crank interval is 4 seconds, and the delay in 
10 implementing the MIC slew is 1 second ( p = 1 ). Assume that the time bias estimated after 

the 64th MIC interrupt is 1 .4 seconds, that there is no frequency error, and that the time bias 
is perfectly estimated. This scenario is illustrated in Figure 15 where the data is identified 
with o*s (the V data): 



15 immediately to 0.4 seconds. This corresponds to point M64+. Note that the time bias is the 
difference between the "o" data and the slope 1 line. Following this, there is a period of one 

second in which the default MIC interval count [n d ) is used. The MIC slew control then 

removes the remaining 0.4 seconds of time bias by calculating a MIC interval count profile 
from equations (28) and (31) . (The quantization effects are ignored in this example since 



After the rounded time bias is removed using equation (15), the time bias drops 
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they will not show up on a plot.) Since the time bias is positive, the unquantized MIC 
interval count is greater than its default value and the true time between MIC interrupts is 
lengthened so that the minimum length constraint is not active. 

Next consider the case where the time bias estimated after the 64th MIC is -0.4 
5 seconds. The true time between MIC interrupts resulting from using equation (3 1 ) would be 
lengthened by 0.4/3 seconds. Therefore the minimum length constraint becomes active and 
the true time in the Kalman crank is increased by one second. The MIC is then aligned with 
the UTC time tics in the interval between MIC interrupts 65 and 66. During this alignment, 
the true time interval between MIC interrupts is increased to 1 .6 seconds. Note that in this 
10 example also, there is a period after the 64th MIC in which the default MIC interval count is 
used.. At the next Kalman crank time (the 68th MIC interval) the MIC interrupt is aligned 
with the UTC tic, but the time bias is now +1 second. Since equation (34) was used to 
propagate the filter estimate of the time bias, the Kalman filter "knows" this. When equation 
(15) is again invoked, the time bias is reset to zero (point M68+). 

15 

Dual-Processor Subsystems 

The GNS software executes on two processors, a tracking processor and a navigation 
processor. In the preferred embodiment, two Mongoose V, Synova, Inc. processors are 
used. They run on an adaptation of the 32-bit MIPS processor architecture. Like the GTA, 

20 the processors are implemented on a rad-hard HR2300 series gate array from Honeywell. 
Figure 16 is a functional block diagram of the dual processor and Figure 17 is a simplified 
layout of the hardware. Both processor hardware configurations are nearly identical with the 
exception that the navigation processor has a command and telemetry interface with the 
spacecraft and the tracking processor has a control and data interface to the GTA. In 

25 addition, the tracking processor has no nonvolatile memory and thus is slaved to the 
navigation processor, which is responsible for booting it up after a reset and continually 
monitoring its health. 
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The navigation processor communicates with the spacecraft via a memory-mapped, 
8-Kbyte, dual-port RAM (DPRAM), which interfaces to PCI hardware on the spacecraft side 
of the DPRAM. The "sharing'- of the DPI AN between the GNS and the spacecraft is 
arbitrated by the GSE by use of interrupts and defined access timing constraints. The 
5 navigation and tracking processors communicate with each other via another 

memory-mapped 8-Kbyte DPRAM and, like wise, arbitrate their access with the use of 
interrupts and defined timing constraints. While the navigation processor has read/write 
access to this entire DPRAM, the tracking processor has read/write access to the upper 4 
Kbytes but only read access to the lower 4 Kbytes. This "read-only" memory provides a 
10 "pseudo-non-volatile memory" for the tracking processor, which is utilized during boot-up. 
The tracking processor communicates with the GTA Via a custom address and data bus used 
for sending control and configuration information and receiving data and status. The 
software executes an extended Kalman filter which is presented by way of a mathematical 
model immediately after the following functional descriptions of the software and its tasks. 

15 

GNS Software 

The GNS software is presented as a functional block diagram in Figure 16. It is 
partitioned between the tracking processor and navigation processor based on functional 
requirements, system topology, and required and available processor throughput. The 

20 navigation processor utilizes an embedded Operating System (OS) such as the Nucleus Plus 
produced by Accelerated Technologies, but only in a limited manner i.e., limited to task 
switching, interrupt handling, and minimal use of events and semaphores. The tracking 
processor does not make use of an OS due to unacceptable interrupt latencies on its high-rate 
interrupt (-667 sec). All navigation and tracking processor software is written in ANSI C, 

25 except for a minimal amount of modularized assembly code. 

Tracking Processor Software 

The major tasks performed by the GNS tracking software are listed in Table 7. This 
software, in conjunction with the GTA, can acquire and track up to 12 space vehicles (SVs). 
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The software supports aidcd-scarch and sky-search acquisition modes for acquiring satellites 
on a channel-by-channei basis. In aided-search, aids generated by the navigation processor 
are used to determine the SV pseudo-random noise (PRN) and the Doppler range searched 
out for potential acquisition. In sky-search, no aids are provided and the system searches a 
5 wide Doppler range for SV PRNs from a predetermined set. For each signal acquired, the 
software uses a third-order, phase-locked loop to track the carrier, and an aided first-order, 
delay-locked loop to track the CA code. A signal monitor is used to determine when a 
channel drops lock. When this occurs, the processor autonomously tries to reacquire the 
satellite with a short search. If this fails, the channel switches to the sky-search acquisition 

10 mode or if aided by the NP, to the aided acquisition mode. 

The software makes carrier and pseudo-range measurements every second for each 
satellite in track, and outputs these data to the navigation processor. In addition, the 
software builds GNS message subframes, validates them using the message parity bits, and 
then out puts validated subframes to the navigation processor. Interrupts generated by the 

15 GTA control the sampling and processing of the data from the 2 tracking channels on the 
GTA. A measurement signal generated by the GTA is used to control processing and 
handling of the pseudo-range and carrier phase data, and subsequent output to the navigation 
processor. 

20 Table 7. Major GNS Tracking Software Tasks 

• Accept and apply acquisition aids 

• Execute sky-search algorithm in absence of aids 

• Control the GTA 

• Track up to 12 SVs simultaneously 

25 • Generate validated GPS message subframes 

• Determine SV transmit time for latched data 

• Execute commands from the navigation processor 

Navigation Processor Software 

47 



BNSDOCID: <WO 0063646A 1_l_> 



WO 00/63646 



PCT/US00/10197 



Table 8 details the major tasks performed by the navigation processor. The core of 
the navigation processor software is the Kalman filter task, which consists of two main parts: 
the Extended Kalman Filler (EKF) and the Short-Term Propagator (STP). The Kalman filler 
updates the state vector and covariance at 1 80-second intervals and the STP generates all of 
5 the data products and outputs them at 1 -second intervals. 

Table 8. Major GNS Navigation Software Tasks 

• Process 1 CCSDS telecommand packet per second 

• Output 6 CCSDS telemetry packets per second 

10 • Output unpacketized data to guidance and control instruments 
every second 

• Generate a pseudo-range table 

• Build OPS message data stores 

• Execute Kalman filter crank every 1 80 seconds 

15 • Output a PVTSE data set once per second; (PVTSE = Position, 
velocity, time, Sun vector events) 

• Generate acquisition aids 

• Generate control for steering 1 PPS signal 

• Generate and output long-term propagation data 

20 The GNS uses a Fault Detection and Exclusion (FOB) RAIM technique. The 

software excludes measurements from the GSE whose residuals exceed a defined threshold 
via the Kalman filter. The GNS uses a current-state EKF to obtain the current best estimate 
of the position and velocity of the TIMED, satellite center of mass. The EKF is essentially a 
simulation corrected by periodic measurement updates. The updates normally occur every 

25 180 seconds and process GNS range and phase measurements for up to 12 satellites in track. 
The state propagation in the EKF contains a Jacchia upper atmospheric density model and a 
gravity model using degree and order 15 spherical harmonics from the EGM96 database. 

The EKF uses a state space containing 45 states, defined in Table 9, to model the 
measurements. The drag error state normally represents the error in the drag model, which 
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is dominated by the unpredictability of the upper atmospheric density. If the parameters for 
the density model are unavailable, the drag error state can be used lo model the whole drag 
effect. 

Table 9. The GNS EKF State Space 



5 Parameter Number of States: 

Position 3 

Velocity 3 

Drag Parameter Error 1 

Clock Error 1 

10 Frequency Error 1 

Selective Availability Markov Process 24 

Integrated Carrier Phase Bias 1 2 

Total 45 



15 The prediction capability of the EKF is primarily determined by how rapidly the 

atmospheric density fluctuates in reaction to solar activity, which will be at a maximum 
during the TIMED mission. 

The clock and frequency error states are used within a control system that drives the 
STA's measurement clock and 1 PPS output to be coincident with UTC epochs. This 

20 latching signal is also used to define the time of validity of the output data products. 

A two-state damped harmonic oscillator driven by process noise models the selective 
availability errors for each S V link. The integrated carrier phase bias is used to model the 
integer ambiguity associated with the phase tracking process. These selective availability 
(SA) Markov and carrier phase biases are reinitialized when new satellites are brought into 

25 track or during reacquisition after a tracking loss of lock. 

After each measurement update, the Kalman filter state is propagated for 180 seconds to 
the next (anticipated) update time. After this, the STP propagates the state for two more 
180-second intervals. These additional propagations are used only for generating real time 
GNS output data products for the spacecraft (position, velocity, time, Sun vector, event 
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notification flags, and data validity flags) and do nol affecl the Kalman filter solution. A 5th 

Hermite position and three velocity points 

degree/interpolating polynomial is fit to the three points/forming the boundaries of the two 

additional propagations. The. interpolating polynomial matches position and velocity at the 

boundary points and is used to interpolate position and velocity to one second intervals. The 

5 interpolated position and velocity are used to generate the five event notification flags: 1) 

primary ground station contact, 2) backup ground station contact, 3) SAA encounter, 4) 

polar region and 5) day/night at sub-satellite point (used for signaling a terminator crossing). 

Approximately every 12 hours, a long-term propagation (LTP) is executed in order to 

predict the primary and backup ground station contacts, and the SAA encounters for the next 

10 60 hours. The duration of the propagation is limited by the required accuracy of the data, 
which is largely dependent on the uncertainty of the atmospheric density at solar max. 
During the LTP, a state vector corresponding to each predicted ground station contact is 
saved and used to generate an orbital element set for the respective contact. All of the LTP 
data products are used solely on the ground (i.e., not utilized onboard) by the MOC and the 

15 instrument Payload Operations Centers (POC). 

Extended Kalman Filter 

The TIMED system thus far described is realized through the application of an 

20 extended Kalman filter. The filter distills information accumulated from received GPS 

signals into current state estimates. This is accomplished by propagation and measurement 
equations which involve nonlinear transformations of the state. A mathematical model 
defining the processes of the program implementation of the extended Kalman filter is 
presented immediately after the following graphic description of the filter. 

25 The time update of the Kalman filter state contains models for the earth's gravity 

field and for drag. However, since the drag model is uncertain, an error state is 
included in the state space. Small errors induced by these and other 
sources are eliminated by a linearization process wherein small errors in 
the state are estimated and added to the whole value state. 
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The filler processes both pseudo-range and carrier phase measurements. The pseudo- 
range measurements offer a noisy (3-5 m with CA code) measurement of slant range from 
the TIMED satellite to each GPS satellite. The carrier phase measurements provide precise 
measurements of line-of-sight velocity which is fully exploited only when the selective 
5 availability mode is off. 

The Kalman filter state space is divided into two partitions which are dynamically 
uncoupled. The first partition includes the states which directly determine the satellite 
position and time, i.e., position, velocity, drag parameter error, receiver clock and frequency 
bias. These states evolve according to nonlinear dynamics. The influence of the receiver 
10 position on the measurements is also a nonlinear transformation. 

The filter estimates current position and velocity in the conventional inertial system 
(J2000). The position and velocity change as a result of gravitational acceleration (as 
provided by a 15 x 15 spherical harmonic model) and the drag acceleration (provided by the 
U.S. Space Command high altitude atmosphere model). Any errors in the drag model are 
15 absorbed by the drag error state. 

The effects of the imperfect receiver clock are modeled by the time bias and clock 
frequency error which have been defined in the Receiver Timing section of this patent 
specification. The time rate of change of the receiver clock is the receiver frequency error. 

The second partition of the state space contains states which do not directly affect 
20 TIMED position and time. They are used to model errors in the GPS measurements. These 
states evolve according to linear dynamics and their influence on the measurements is also 
linear. 

These selective availability states are a set of second-order Markov processes, one 
for each link, which model the effect of selective availability on the pseudo- range and 
25 carrier phase measurements. A second order process is used with an undamped natural 

frequency of 300 seconds, damping ratio of 0.4, and a steady state standard deviation of 75 
m. When a GPS track is initiated, the selective availability state standard deviation is 
initialized with the steady state standard deviation. 
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The carrier phase hias states model the unknown cycle ambiguity in the phase- 
derived range. The carrier phase bias state standard deviation is initialized at Le9 when a 
new track is initiated, or when loss of phase-lock is detected. 

The Kaiman filter can be initialized in one of three ways: 1) The normal mode of 

5 initialization through pre-loaded orbital insertion data. The separation position and velocity 

in an ECEF coordinate system are loaded by command prior to launch. When the separation 

initial conditions 
event is detected, the Kaiman filter transforms the /to J2000, initialize its state and 

covariance, resets the crank interval counter, and begins cycling. 2) through a GPS doublet. 

In this mode, the GPS point solutions are used to establish the initial position, velocity and 

10 associated covariance. This process assumes that the point solutions are close enough in 
time that the TIMED velocity can be regarded as constant . The Kaiman filter 

uses the GPS doublet to reinitialize this state if it detects that its state vector is inconsistent 
with a large number of GPS measurements. 3) through a commanded upload. In this mode a 
NORAD tAvo line element set (compatible with the SGP4 orbit propagator) is transformed to 

15 a J2000 state vector and covariance by the ground system. This mode is intended for 
unusual cases in which the GPS cannot function. 

The filter time update propagates the state and covariance between measurements. 
The filter time update is done immediately after the measurement update so that the MIC 
corrections can be spread evenly across a Kaiman crank interval. This requires that the 

20 Kaiman filter measurement updates be evenly spaced. 

Every 180 seconds, difference measurements are formed on the basis tracker and 
GTA data latched at the corresponding MIC. The difference measurements are the pseudo- 
range and carrier phase minus their predicted values. At this time, measurement sensitivity 
is calculated. The filter state and covariance are updated by the measurement update 

25 process. The Kaiman filter residuals are monitored to verify that there are no huge 

discrepancies between the Kaiman filter state and the measurements. If there are, the filter is 
reinitialized using the GPS doublet. After the measurement update, the GNS output 
products are generated and a check is made to determine whether the long-term propagation 
should begin. 
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The independent variable in this state propagation is receiver time. An increment of 
receiver time differs from an increment of true time because of the receiver oscillator 
frequency error. The position, velocity, drag, clock and fractional frequency error states (the 
dynamic states) are coupled dynamically to each other by a system of nonlinear ordinary 
5 differential equations, the derivations of which commences with equation (35). A transition 
matrix that propagates small errors in these states is generated using the gradient of the 
dynamic state derivative with respect to the dynamic state. This discrete time process noise 
matrix for the dynamic states is also found by integrating its o.d.e. which involves the same 
gradient. 

10 The system of nonlinear ordinary differential equations describing the state dynamics 

is solved by using an 8 th order Runge-Kutta Fehlberg algorithm. The performance of the 
filter begins to degrade because of measurement thinning if the step size is increased beyond 
180 seconds. Because propagation is from measurement to measurement, a fixed-step size 
version of the algorithm is used. The RFK8 algorithm requires 12 function evaluations per 

15 step. The subroutine that implements the algorithm integrates vector (not matrix) 

differential equations. Therefore, the discrete time process noise matrix must be remapped 
into a vector as well as its derivative. 

The selective availability (SA) states and the phase biases (the "sat" states) obey 
linear dynamics. The separation between measurements is constant, allowing the continuous 

20 part of the transition matrix and discrete time process noise to be precomputed. The "sat" 
states are ordered according to tracker channel. The transition matrix (and the discrete time 
process noise matrix for the "sat" states) model the effects of a new GPS satellite in a given 
tracking channel. 

A full state ("dyn" + "sat") transition matrix is formed from the block-diagonal 
25 concatenation of the "dyn" and "sat" transition matrices. The "dyn" partition is computed 
by solving an o.d.e. containing the "dyn" state derivative gradient. The covariance factor is 
updated for deterministic effects by premultiplying it by the full state transition matrix using 
a sparse matrix multiply routine. To accomplish the stochastic update of the covariance 
factor, a vertical concatenation of the deterministically updated covariance factor and the 
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full-state discrete time process noise factor is formed. This array is then decomposed into an 
orthogonal matrix and upper triangular matrix using a sparse QR decomposition algorithm. 
The updated covariance factor (transpose) is found in the upper part of the upper triangular 
matrix. 

5 The raw pseudo-range and carrier phase measurements are formed from the 

following quantities: 1) fractional chip count (pass-through from GTA); 2) chips into epoch 

(pass-through from GTA); 3) epochs into bit (from tracking SW bit sync process); 4) high 

order bits of Z-count or TOW (Time Of Week) (from message); 5) GPS week (from 

message); 6) whole cycle count; 7) fractional cycle count; and 8) receiver time. 
10 Predicted pseudo-range and carrier phase for each satellite in track is then calculated 

from the receiver position, the SA and phase bias states of the GPS message ephemeris. 

Since the predicted pseudo-range is a function of the satellite position at the predicted 
transmit time (which in turn depends on the predicted pseudo-range) this process is 

iterative. The Kalman filter processes the difference between the raw measurements and 
15 their predicted values. The measurement sensitivity depends on the unit vectors from the 

TIMED satellite to the GPS satellites in track. They are produced in the calculation of 

predicted pseudo-range. 

In the Levy square root covariance Kalman filter derivations, the array on the left of 

equation (67) is formed from the measurement noise factor, the time-updated covariance 
20 factor, and the measurement. This array is factored into an orthogonal matrix and an upper 

triangular matrix containing the measurement update covariance factor, and the W matrix 

which is used to compute the Kalman gain. 

After every measurement update, the filter launches a low precision state propagation 

to fill a table containing position, velocity, sun-vector, event notifications, and time for the 
25 next 180 + seconds. In this mode the orbit propagator may only use the J2 part of the 

gravity module to improve execution speed.. The event notifications are: 1) entry into the 

South Atlantic Anomaly (SAA) or Polar regions, 2) ground contacts and 3) terminator 

crossings. 
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Every 12 hours, the Kalman filler launches ;i high precision orbil propagation usiny 
its last computed state vector. This trajectory is monitored for contact events and any SAA 
entry. When a contact event is found, the trajectory points near the AOS are used to define 
an element set appropriate for used in the NORAJD SGP4 algorithm. After every 
5 measurement update, the Kalman filter produces acquisition aides for the tracking software. 
The aids contain Doppler and Doppler rate for all GPS satellites based on the current 
TIMED state vector and current almanac ephemeris. 

To avoid missing an event because of the granularity of the propagated orbit, the 
event detection and prediction is based on polynomials fit to functions of the propagated 

10 points. For the contact events, a polynomial is fit to the elevation of the TIMED satellite as 
seen by the various ground stations evaluated at the trajectory points. The roots of this 
polynomial are the AOS (going +) and LOS (going -) times. If AOS - LOS>5 minutes, it is 
a valid contact. For the South Atlantic Anomaly, a polynomial is fit to the latitude and 
longitude polynomials of the sub-satellite point. The points at which the latitude and 

15 longitude polynomials cross the latitudes and longitudes of the SAA region is calculated. If 
the satellite is in the SAA latitude band at the same time it is in the longitude band, an SAA 
event is flagged. The polar region is done in a similar fashion, except only the latitude is 
used. For the terminator crossings, a curved is fit to the inner product of the sun-vector and 
the unit normal to the reference episode at the sub-satellite point. The roots of this 

20 polynomial define the terminator crossings. 

Kalman Filter Implementation 

The following presentation mathematically describes the processes carried out by the 
software in implementing the various routines used to accomplish the tasks of the extended 
25 Kalman filter which forms the nucleus of the GNS. 

The state space for the TIMED GNS Kalman filter is: 
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:(/) = 



"w(r(/)) 
a(t) 
T(t) 
sit) 
sa(t) 
bit) 



9 + 



(35) 



where w(r(t)) is the position and velocity of the TIMED center of mass (CM) in the 
conventional inertia! system (CIS): 



w< 



<r(/)) = 



r(rit)) 
v(r(f)) 



(TIMED inertial position & velocity). 



6x1 



5 Note that w(r(/)) is defined as a composite function: w(«) is an explicit function of true 
time. The evolution of is governed by Newton's second law, which involves 

derivatives taken according to true time (not receiver time). The units of are meters 

and meters per second. The drag error parameter is defined as the ratio of the true value of 
the parameter group p • C D • A/(2 • M) to its nominal value: 

, . ( p-C D -A ) ( p-C D -AY omina ' 
10 ait) s TT— - ~ (drag parameter). 



2- M ) V 2 • M 



The state a includes the effects of errors in density, drag coefficient, projected area, and 
satellite mass. The first term on the right is the true parameter group, the second is the 
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parameter group predicted by the density model, with the nominal spacecraft ballistic 
coefficient, mass, and an average projected area. 

The clock related, states T(t) and s (t) and are defined according to: 

T{t) = / - r{t) (Clock bias) (36) 



fo e/ -fo(r) 



s{r) = — — -j^f (Fractional frequency bias). (37) 



The effect of selective availability is modeled as a two-state Markov process for each 
satellite. The vector sa{t) contains two states for each satellite that can be tracked for a 

total of 2>7 max states: 

sa(t)=\sa x {t) sa 2 (t) • • • sa 2 ^ x {t) ■ (38) 

10 The odd elements of sa{t) represent the effect of SA in meters on the pseudo-range and 
integrated carrier measurements. The even elements are the (receiver) time derivatives of 
the odd elements. The vector b{t) is an array of phase derived range biases to account for 
the arbitrary initialization of the integrated carrier phase measurement: 

b{t) = \bAt) bAt) • • • b„ ,(r) b (t)] (Phase derived range biases). 

I » * max 4 mi\ J 

15 The equations defining the forming and processing of the Kalman filter 

measurements use the index k to denote times at which a Kalman filter measurement is to 
be processed, which is every m th MIC interrupt. Previously, k was used to denote the 
MIC times . 

The pseudo-range measurement begins with the raw data where the GPS signal 
20 coming into the TIMED antenna at true time r is : 

^(r)= c4P/^\r(r))-m(r(r))-sin(2^ v (r)-^ C, (0) (39) 
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where PRN V is the PKN (p.seudo random noise; nunu>er uFihe satellite in tracking channel 
v, CA(PRN V ,/ r ) is the CA (coarse acquisition) code transmission for the satellite as a 

function of the signal-indicated transmit time f v (r), m(/ v (r)) is the data modulation of 

the GPS transmission, sm(2xf LX t v (t) - ^^ Cl (r)) is the GPS carrier, / t| is the 

1575. 42 

5 nominal GPS carrier frequency v Mhz, and pr^f (r) is the phase retardation of the 

GPS carrier (relative to the code) caused by ionospheric dispersion. This is equal to the 
difference in the phase and group delays: 



t iono 



W-^Wr.)-/;^.)) (*» 



where is the group delay in seconds and I v p {? k ) is the phase delay. The GNS 

10 receiver generates a local replica of the GPS CA code, which is synchronized with the 
received code when the receiver has locked-up: 

CA(PRN\t v (r))= CA(PRN\t{r)) 

At receiver time t k = t(j k ^) y the local code phase is latched, this is the measurement 

interrupt clock (MIC) event. The value of the transmit time in seconds since the GPS 
15 standard epoch in terms of the local code phase and quantities in the GPS message is: 
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t x k = /_ int* + / float X j. 
t_ int A v = GFS_week\ • 86400 ■ 7 + 7(9 • 6 

bits_n2_ sframe x epochs_n2_ bit x 



1 .float; = 



+ 



50 1000 
chx_cacode_ phase + chx_codenco_ phase /2 

1023000 



16 



(41) 



Here t v k is the receiver estimate of the transmit time of that point in the CA code that the 

local code was correlating with at the k th MIC. It is composed of an integer part and a 
floating point part which can never exceed 6 seconds. The quantities bits _n2 _ sframe x 

5 and epochs \n2_b it x represent the number of completed bits into the current subframe and 

the number of completed code epochs into the current bit for tracking channel x . They are 
maintained by the tracking software after it has achieved bit and frame synchronization. The 
quantities chx_cacode_ phase and chx_codenco_ phase are read from the registers on 

the TIMED GTA chip containing the full chip count and the fractional chip count for 
10 tracking channel x . The variable GPS^week (the current GPS week) is initially set by 

the contents of word 3 of subframe 1 in the GPS message and is incremented by one (and 

seconds Ji2_ week is decremented by 7 -86400) if seconds^ 

n2_week= T_float x k exceeds 7-86400. TOW is the number of 6 second 

counts from the beginning of the week to the beginning of the next subframe, taken from the 
15 handover word of the last completely assembled subframe. 
The pseudo-range is then constructed as : 
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U \ 1 (42) 

.= c\(t k -t_\nV k )-t_floatl\ 



where c is the speed of light, and t k = MJCjime k + / in1 is the (controlled) receiver time 
corresponding to the MIC. Note that the large integers / A and ?_ int^ are differenced 
before subtracting the floating point number t _ float k so as not to lose precision. The 
5 model for the indicated receive time is : t k - T k + T k . 

The signal indicated transmit time is the transmit time derived from the modulation 
on the GPS signal. The model for the signal-indicated transmit time is: 

/;= r; + r(r;)+-5a:( r ;) (4 3) 

The receiver indicated transmit time is the signal indicated transmit time corrupted with 
10 receiver tracking errors: 

i; - /; + se k 

(44) 

sa:\T;)+6t: 

c 



= r; + r(r;) + ^^;(r;)+^ v 



where t k is the signal-indicated transmit time and St k is the tracking error: z k = T v (t k ) is 

the true transmit time of the signal being tracked in channel v at the k th MIC interrupt, 
2" v (r* ) is the GPS satellite clock bias, sal'^rfy is the selective availability process, and 

15 the receiver tracking error, StJ : , is assumed to be white at the MIC interval. The true 
transmit and receive times are related by: 
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r - r- - f 
/; * J/m /// y (r 



in 



where v (r) is the group speed of the GPS transmission as a function of position, and 

true 

where the integration path (PATH) begins at the GPS satellite position at the/transmit time 
( r^(r A v ) , and ends at the position of the receiver antenna at the receive time 

5 (r on '(f t ) = r(r k ) + T^.(r A . )L 0 „, ) . Here r(r k ) is the TIMED center of mass position i 

the CIS system at the true receive time r k , r" is (r v t ) is the true position of the GPS satellite 
in track by channel v at the transmit time r v k , L a „, is the CM to antenna phase center lever 
arm in the TIMED body system, and T^,{r) is the transformation from body to inertial 

coordinates obtained from the attitude system. 
10 Assuming: 

111 



1 r 



1+ 1- — 

V c J 



the free space range delay and the ionospheric delay are separable, so that the time 
difference may be modeled as: 
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J PATH y (r 



r hAThr 1 ^ ipATH\ /- 11 



C J PATH\ 



= 4'(r-A'(r))-r-(r| t if 



1- 



\dr\ 



= ^(0 + j;(r k )/c 



where A v (r) is the free-space range delay (in seconds) for the v' h link, and 

1 



r(r) l i i v * (r) l 



= -f (l-«(r))kr| 

C ^ PATH " 



(45) 



is the ionospheric range delay (in meters) at the true receive time. Here n(r) is the index of 
5 refraction as a function of position. The free space range delay is given implicitly by: 

1 



A'(r) = -r'(r-A-(r))-(Kr)+1- (r)^) 



1 



= -|r v (r-A v (r))-r°"'(r)| 
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where r ani (r) = r{r) + T™ dy {j% a ,u is the position of the TIMED GPS antenna. The 
model for the raw pseudo-range measurement is: 

Pl = c\t k -t;) 

= c\r k + 7i. - r ; - r(r;)- ^. v (r;)/c- 

/ N (46) 

= c-(A v (rJ+ - r(r;)- + /;(r 4 )+ *,"(r 4 ) 

where a second order markov random process in the range domain keyed on receiver time, 
sa v (r k ) , is used to model the true selective availability process - sct*{rl} (note the sign 
reversal). 

The Kalman filter pseudo-range difference measurement is constructed as the 
difference between the measured pseudo-range and pseudo-range predicted from all 
measurements processed up to the previous time : 

- pI - pl/k-i (47) 

The predicted pseudo-range is given by: 

Pk/t-\ = c-(av. + 7U_, - c„fc v )+ ^ v (v-.)) 

15 where: 
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10 



= estimated free space range delay 
Jg^essitk) s group delay from GPS message 

T k/k _ x = estimated receiver clock bias (controlled) 
^mewO**) ~ GPS satellite clock error derived from message 

= a l + ^. v (siw(f;)- r 0c ) + ^( S iw(f;)- ; 0c ) 2 

siw(f; ) = seconds into GPS week corresponding to t k 

a] = Clock compensation the params for sat v (all from current message page) 
t 0c = epoch for clock compensation (seconds into week) 

estimated selective availability markov (range) 
= salj~^ (See equation (38)) 



The estimated range delay is defined implicitly by: 

c 



(48) 



where: h r k , k _ x + T^(,, )L anl 

T k i k _ x = estimated true receiver time 



15 =h-T k/k _ x 




BNSDOCID: <WO 0063646A*1_L> 



WO 00/63646 



PCT/US00/10197 



f A y"_, is the estimate of TIMED GPS antenna position in the CIS coordinate system at the 



receive time and r v [f k j kmm} - A^_, j is the position of the GPS satellite tracked by channel 

v at the estimated transmit time. (The procedure used to solve equation (48) for 
follows.) 

An iteration is used to find the solution to the light-time equation: 

o=#(av.) 



where: 



<f(A)= A-^r*(f^. I -A)-f t/4 . 1 

C 

r v (r) = IT (r)r"(r) 



r c v w (r) is the position of the v' h satellite in the earth-fixed CTS frame calculated from the 

10 GPS message ephemeris. 

Using an initial guess for A : 

1 



A(0) SS -r"(f 4/4 . I )-f it/4 . 1 



(49) 



the iteration proceeds as: 



(50) 



15 until jA + l) - A ( y )| < tol . Typically, tol - l.e - 1 1 seconds. At convergence set 
A*/*-« = A 0')- 
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In the Kalman filler, the pseudo-range difference measurement is determined by 
substituting into equation (47): 

y;(k) = c(A v {r k ) + T k - r(r;)- st:)+ /;(rj+ 
= c(^ + ^) + ^ 2v -' + ^ 



this assumes that (*; )= j;(r k ) and C w fe)= r (^) (no errors in the 

5 ionospheric compensation or GPS clock compensation in the GPS message data). 

$Pk ~ c ° £ { k * s a w hite noise process representing all high frequency delay-locked loop 

errors, and SA v k = A v (r) - A v ^_, is the error in the estimated free space range delay. The 

solution for SA\ in terms of the Kalman filter error state follows: 

The errors in position and velocity are synchronous, therefor: 

V-, = r(r(/,))-^r 

/ \ (51) 

Assuming that the attitude system provides an orientation matrix that is a perfect real-time 
reference and the lever arm is perfect: f™ v (/, )L anl = T^.(r, )L 0 „, so: 

= r""{r k )-Sr k 
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The true quantity minus the error is substituted for the estimated quantities in 
equation (48): 



Where: 



0e A*(r,)- 8A\ - ±\r*(t k - f k/lc _ x - £S{r k )+ 8A\)- r<"*(r,) + 8r k \ 
= tr{r k )-SL\ - \\r v {t k - (T k - 5T k )- A"(r,) + 8L\)- r~{r k )+ 8r k \ 
= A v (r k )- 8A\ - ±|r'(r 4 + ST k - A v (r,) + 8A\)- r-(r 4 ) + <?r,| 
^ A'(r 4 )- - i|r"(r 4 - A v (r)) + ^(r, - A v (r,)) (^. + 8A\)- r*"'(r,) + ^ 



Where: 



e t = 



.,)- 



~ant 
T k/k-\ 



V 



r 



(from GPS ephemeris). 



But A *(r)-4"(r t - A v (r t ))-r""(rj| = 0 , so 
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Solving for ^A, yields: 



SA\ = 



{& k )'\* k + V-8T k ] 

(e;)V 



' dA v ~ 


6T k + 


~ dA v ~ 


dT 


k 





(52) 



where: 



<?A V 
dT 



dA v 
dr 



fey 

+ fe)'v v 



(53) 



(54) 



The measurement model for the pseudo-range observable is then 



y v M)= c - 



dt£_ 
dT 



+ 1 



\ST k + 



dA v 



k j 



+ tea]"-* + dpi 



(55) 



leading to the concatenated measurement model: 



68 



BNSOOCID: <WO 0063646A1J_> 



WO 00/63646 



PCT/US00/10197 



y. 



,(*) = 



h m . o H r o n sa o 



Sw 

5a 
ST 
8s 

Sh 



(56) 



where 5p k is modeled as white Gaussian noise with (dp ] = dia^0 2 pk ) and where the 



measurement sensitivity partitions are: 



'dr 



or 



0 
0 



db n " 



0 



(-57) 



[*'/*] 



as defined in equation (54) 



69 



BNSDOCID: <WO 0063646A1J_> 



WO 00/63646 



PCT/US00/10197 



+ 1 



+ 1 



H r = c- 



(58) 



+ 1 



and 



as defined in equation (53) 



H so = [l 0 1 0 • • 1 o[ 



(59) 



The raw phase measurement in cycles is constructed as : 

chx_ cyc_ cnt lower 



c/ix_ cyc_ upper + 



.16 



(60) 



where dix_cyc_cntjupper, and chxjcycjontjupper are the registers on the TIMED GTA 
chip containing the whole cycle count and the fractional cycle count for tracking channel x. 

10 In deriving the truth model for the raw phase measurement, the oscillator phase as a 

function of true time is given by: j^(r)= lLnf™ f {x+ T*(t))+ $ set , where T* is 
the uncontrolled time bias. 
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The GPS signal coming into the antenna is : 

s„(t)= c4P/^\/ v (r))-m(^(r))sin(2^ 1 ^'(r)-^ / 0 c, (r)) 
If the phase of the signal after code and data wipe-off is ^ R ( r ) : 

5 as developed by equation (43), t v k = r A v + T v [t^ + — sa*{rl^ and 

r,-r; = A v (r,) + ^/;(r,) 



then: 



KM- 2^,(r; + r(r;) + - W 



= 2^ 



The ionospheric phase delay relative to a wave propagating in free space is : 

io /;(*-) = -/;(r) 

The minus sign indicates the wavefront arrives before a wavefront propagating in free space. 
Referring to the definition P/^ CI (r) in equation (40), the phase advance relative to the code 



is: 
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where /^(r^) the ionospheric group delay defined above. Thus: 

^0^Jr. + r(r.)-V(r.) + (^) +JB ' W j/ 



The downconverter-sampler effectively mixes the received signal with a 
5 k x multiplication of the fundamental oscillator. The combined effect of the downcon version 



and sampling process result in a value for k x of: 



k. = 140 



( i n 

, 1 + — + — 
V 10 AS) 



4 

+ — 
7 



downconvener eflecl sampling effect 



The continuous-time representation of the output of the downconverter-sampler is: 
s IF {r)= LFP{ 

- ^cos[^(r)- k } y(r)] 



10 where LFP{x(t)} represents a low-pass fi herring of the process x{t) . When locked, 



the 



phase tracking loop will match the total phase of this signal, modulo 2n . The model for the 
raw phase measurement (in cycles) is then: 
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[<Pk(t) - k x \j/ (r)] 




2/7/+ A /v " + <ty r 



where jV v is the cycle integer ambiguity which is constant as long as channel v retains 

phase lock, and 8<j) v represents the phase tracking error, which is regarded as a white noise 

process at the Kalman crank iterval. The GNS tracking processor samples this process at the 
5 MIC interrupts and performs the following deterministic transformation on the raw phase 
samples to provide the Doppler sense that is conventional to the GPS post-processing 
community (i.e. consistent with the RINEX file definition), and to eliminate the large 
variation in phase caused by the known part of the zero Doppler offset: 




(61) 



zero Dopplcr offset 



10 



This transformed phase data is what crosses the boundary between the tracking 
processor and the navigation processor. The model for the phase data is derived by 
combining the above results: 
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and A u = c/ f u is the GPS LI wavelength. From equation (4), t*(r k )= r k + T m (r k ) 
leading to: 
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^(r,)=/ 11 (r(r,)^r(r t ) + A'(r,)) + ("' ( '' ,_/;(r % + 6*-* ) »* 

= /„(*,)- u(r,)- r (r,) + A'(r.)) + ^"^'^ '^'^ + 6' - 

when: 

7tO-r(r t )=r-rM-r-rM 

= (63) 

and where the biases have been absorbed into b x 

The observable processed by the Kalman filter is the difference between the phase 
measurement and its predicted value based on the current state estimates: 
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The estimate of the phase measurement is given by: 




The phase difference measurement is given by: 

(64) 



5 
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Because of the similarity between equations 64 and 55, the concatenated 
measurement y„(k)= J^UXj^U) - ■ * ') ; ^T (*)] can be in the same form as equati 



(56): 



h m o n T o h m o 



5xv 

5a 

5T 

5s 

Jsa 

5b 



(65) 



where 5<j) k is modeled as white Gaussian noise with cov{5(f) k ) = diag[cT ^ k j 
The concatenated measurement matrices is defined as: 

V 

H 



y* = 



H 



(66) 



cov 



{Sp k ) o 
0 cov{5</> k ) 
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The measurement update proceeds as follows: 
R k = R k /2 R' k /2 is factored where R^ /2 is the upper triangular. Note that since R k is 

diagonal, Rf = diag^R 

Q - R decomposition is preformed 





*R' /2 

*^k 


0 




"Rf 




5 


p'/2 
k/k-\ 


H l P' /2 


= T 


0 


r k/k_ 



(67) 



10 



15 



where P^^ r is the factor of the prior covariance, B^ 2 iS f act or of the residual 
covriance, and = W Jt B*^ 2 is the Kalman gain. The posterior covariance factor is 

k/k ' 



The measurement update of the state is then: 

<*x k/ * = K^y, (68) 
Since: 
&k/k-i = 0 

The whole-value dynamic state is re-initialized as: 
X k/k = X k/k-\ + ^**/* 

The evolution of the position and velocity states depends on the drag and clock 
states, but not the SA states or the PDR range biases. The state dynamics can therefore be 
decoupled by partitioning the state as follows. 
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w(r(/)) 

a{t) 

At) 

AO 

sa(r)" 
b(t) 



The time update is then represented as : 



x 



sat 



.('...) 



^( X dyn(^)^*^* + .) 

^(t k+u t k )x sat (t k \ 



(70) 



The update of X sat can be written in terms of the transition matrix O l(t k+x ,t k ) 
5 since its dynamics are linear. The vector ^(x dyn (^t k ,t k+i j is the solution of : 



a?/ 



x dyn (^)= f dyn (x dyn (/),/) 



integrated from t k to + 1 with x d (/^ ) as the initial condition. 

By using the chain rule, the position and velocity partition \v(r(f)) is evolved 
10 according to: 
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d\\ _ dw(r(/)) dr 
dt " dr dt 

= f w (w(r(/)) f £r f r(/))r(r) 



(7!) 



Where r(f) is evaluated as described with respect to equation (11), and where the 
rate of change of w(#) with respect to true time is f w (w(r(/)),ar , r{t)) . The term r(f) 

is the ratio of a true time differential to a receiver time differential and arises because the 
5 system is integrating with respect to receiver time. The explicit relation for the true time 

derivative f w (w(r(/)),tf, r(t)) is: 



ar 



r(r) 
v(r) 



Mr) 
dr 

dvjr) 
dr 



= v(r) 

= T c ^(r)-g cw (r cu (r))- 



(72) 
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where r„( r) = T£( r) • r(r) is the CTS position of the TIMED CM, g CM (r m (r)) is the 



acceleration of gravity, T™(r) is the transformation from the conventional inertial system 

to the conventional terrestrial system (see "GPS Theory and Practice" by B. Hoffrnann- 
Wellenhof, et. al., published by Springer-Verlag 1992, the contents thereof which is 

5 incorporated herein by reference), p(r( r), r) is the mass density from the atmospheric 
model, and fi 0 is the nominal TIMED ballistic coefficient, defined as ■ A . The 

earth-relative velocity coordinatized in the CIS frame is given by: 

v re/ (r)= v(r)-a> e x r(r) 

where CO e is the angular velocity of the Earth in the CIS frame. Note that the explicit 
10 dependence on receiver time is suppressed with the understanding that r= r{t) . The time 
rate of change of the clock error is given by the following form of equation (12): 

u(t)- s(t)_ 

(73) 



f(t) = 



1+ u{t)- £{t) 

= {u{ t )-s{t))m 



The time rate of change of x dyn is then: 

■f w (w(r(/)),a,r(/))-rW" 
0 

(u{t)-s{t))r{t) 
0 

1 



f dyn (x dy n(0^) = 



fit) = 



\+u(t)-s(t) 
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The PDR biases are not dynamically coupled with the SA states, therefor the 
transition matrix <t>2^t k+l i k ^j can be partitioned 



as: 



^ sal + "" 



(75) 



The PDR bias states are true biases, they have no dynamics. However the transition 
5 matrix must take account of rising and setting satellites. (Loss of lock is modeled by 
increasing the covariance on the PDR biases to infinity.). To do this, a list of the PRN 

numbers is constructed for the satellites in each tracker. At receiver time t, PRN{t) i is 

the PRN number of the satellite being tracked by receiver channel / , for 

/ G {l,2,...,/7 max | . If channel / is not tracking any satellite PRN(t) t = 0. The selector 

10 matrix Sl(/ ; . +l5 / A ) is defined as: 

Sl(w*) lt/ = s(pRN(t k J r prn(t k ))} 



(76) 



Where S(ij) = < 



1 /'=./> 0,y*0 



0 otherwise 



ij e {l,2,3...,/7 max } 



For instance, assume a scenario with 3 tracking channels (« max = 3) and PRN 10, 

14, and 22 are being tracked by channels 1 , 2, and 3 at t k . At measurement time i k+i , PRN 
15 14, 10, and 20 are being tracked by channels 1, 2, and 3: 
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PRN(t k ) = 



10 14 22 



14 10 20 



si(w*)= 



0 1 0 

1 0 0 
0 0 0 



Sl(/ < + 1 ,t k ) connects the satellite identities at time t k to the identities at time 



t k+i . Then: 



(77) 



For each satellite, SA is modeled as a second order Markov process: 
The solution can be written in terms of the transition matrix as: 



d_ 
dt 



n 2 



1 0 



5,W 



^-s(/)= F mih .-s(/) 
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Written in terms of the transition matrix, the foregoing is: 
4J = ^h + i-';)'4) 
where i. 0 is given by: 

> h ) = a 0 I = ff , • F mkv ■ A / (2x2) 

(78) 

or 0 = -e u (w- sin*y - # • cosco)/q) 
a x = e" • sino) /a) 



5 The full SA state vector is the concatenation of « max copies of the per-link SA states, so that 
the full SA transition matrix is block diagonal with 0 miv [t k+i , t k ) on the diagonals: 

*■*,('*♦,.'*) = diag{<} mkv {t k ^t k )] (2-n max x 2-/7 max ) 
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where t~ k ^ is a fictional time prior to updating the satellite states to account for the addition 
and deletion of satellites. Bookkeeping is accomplished with the selector matrix defined in 
equation (76) by segregating the S A states into terms associated with range (sa_l) and 



those associated with range rate (sa_2) : 



sa 
sa 



_ 2 (^.). 



= P sa 



where P is a (constant) permutation matrix given by: 



P = 



E_ = 



e = 



E, E 2 • E m • • E r 



e m 0 

m 



0 



0 0 • • 0 1 0 • «0 

m lh component 



(2« m ax X2 w m ax) 

(2 • « mav x 2) 

\ max / 



(80) 



10 



The addition and deletion of satellites is accomplished with the selector matrix 
defined above: 

sa_l(Ci)l [S1(W*) 0 
sa_2(C,)J 0 Sl(r, HlJ /,)J. 



sa 
sa 
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P is a permutation matrix, there for its inverse is its transpose so the SA states arc 
returned to their original order by premuitiplying by the transpose: 

"sa_l(C,)" 
.*a_2(/; +] ) 



sa('* + ,)=P' 



The complete update to the SA states is: 
5 sa(/, +1 )= O sa (/, +I ,/,)-sa(/,) 

where: 

3>sa , h ) = P' • S2(/, +1 , t k ) 

si(w*) o 



o si(w 4 )_ 



(81) 



The time update of the extended Kalman filter state estimates is exactly the same as 
the evolution of the Kalman filter state described in the foregoing truth model, i.e.: 



10 



satijk + l/k ) 



^ sat (h+ 1 > h ) ' *sat 



86 



WSDOCID: <WO_0063646A1_I_> 



1 



WO 00/63646 



PCT/US00/10197 



The TIMED Kalman filler is a square root covariance implementation wherein the 
state error covariance P = ccv( Sx) is maintained in factored form as P = P l/ " ■ P l/2 

where P' /2 is upper triangular and P 1/2 = (P' /2 )' - The the covariance factor update is 
achieved by postmultiplying the transpose of the transition matrix: 



(82) 



where Pj/^ /k is the state error covariance after the deterministic update and prior to the 

stochastic update. 

The full transition matrix used in equation (82) is: 

0 O ».('* + !.'*) 



10 where O sa( (t k+l , t k ) is as computed above, and ® dyn (t k+x , t k ) is the solution to: 



(83) 



between £ = t k and£=r t+) . The initial condition is <D dyn (t k , t k ) = I. F dvn is the 



gradient of f d with respect to x dyn = [w'.flj.f]' as derived below: 



87 



BNSDOC1D: <WO 0063646Al_l_> 



WO 00/63646 



PCT/US00/1O197 



<7W 


da 


dT 


ds 


0 


0 


0 


0 


0 


0 




df 


0 








de 


0 


0 


0 


0 



The partitions of F dyn are: 





d\\ 








da 


*(Vr) 






dT 







Si 



= 0 



r = - 



<?r dT dr 
= i 
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In the last two rows above: 



dr 



ds de (l + u - z) 



(l + u - z) 



- T 



The above results are substituted into equation (84): 



F 

A dyn 





di„ . 


■ T 






3a T 


dr 




0 


0 


0 


0 


0 


0 


0 


-r 2 


0 


0 


0 


0 



(85) 



3f / 

5 where r = (l + u{t) - f (/))"' and f u . is given by equation (72). The /fa partition 



is: 
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or 


or 




~dr 


ov 


ow 


dv 


ov 




Ik 


Ov. 


or 


0 




dr ~ 




dr 






dv ~ 




ov 
or 


d < 
or 





2-A 



{l+a)\v rel \v 



rel 



ryCIS 
- '■cts 



(86) 



T^r)^T-(r) + 4^0^) 
^cs 2- A) 



v rg/ -v re/ I+ v re/ (v re/ ) f 



re/ 



{^ X } 



<7V 

~dv 



d_ 

(TV 



T c "(r)g (r )_^ r -' r ) 

Vtt V * ' *=>cts\* cts ) 



0 + *)| 



2-A 



2-A, 



v re/ v re/ 



re/ 
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In die above differentiations, variation of density wilh posiiion are noi considered. The 

df / . . 

/da P art,tlon ,s g ,ven py- 



da 



dr 

da 
dv 

I da J 



0 



P 



2 k 



v rel • v reI 



df , 

The /g T partition is: 



dr 



dt 
dv 

dr . 
v(r) 
v(r) 



(87) 



where v(r) is as given in equation (72) and the important terms in v(r) are given by: 
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v(0=^jT~(r)g„,(T»(r)r(r))-4^(l + «)kiv"'' 



2-A 



^ v re/ (v re/ ■ v re/ ) + v fe/ (v re/ • v^')^ 



2-A) 



V 



rel 



J 



" tj(0-S.(r.X(:)^fer t C«) 



/>( r c*> 0 f v re/ (v re/ • v' e/ + v rel (y rel ■ y rel Y 

+ L 



2-A 



rel 



(88) 



where the time rate of change of the Earth relative velocity is: v rei = v - m ex V 



(89) 



The effects of process noise on the covariance factor are accounted for by the 
following Q-R decomposition: 



p</2 

r k + ]/k 

0 



p</2 
A k + \/k 

QD' /2 



(90) 



where T is an orthonormal matrix and QD is a discrete-time process noise covariance 
matrix for the entire state, which is partitioned as: 
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0 


0 


QD 


0 




0 




0 


0 





(91) 



where QD dyn = COv(<£x dyn ) is the error co variance of the x dyn partition of the state. 



QD dyn is the solution of: 



4f 



QD dyn = F dyn (<f) • QD dyn + QD dyn ■ F d ' yn (<f) + Q dyn 



(92) 



10 



5 between £ = t k and q = r &+1 . The initial condition is QD dyn = 0 . The square root 

QD 1/2 is given by: QD 1/2 = US 1/2 where QD = USU' is the SVD factorization of (the 

symmetric) QD . ( QD 1/2 obtained in this way is a symmetric square root, not lower 

triangular. Although lower triangular square roots are baselined in the Levy SRCF 
algorithm, any square root will do) 

The process noise for the dynamic states ^Q dyn ) is essentially a tuning parameter. 

The physics of the problem dictate that: 1) There should be no process noise driving 
position. 2) There should be enough process noise driving the velocity to cover all of the 
unmodeled accelerations (third body effects, Earth tides and so on). 3) The process noise on 
15 the drag parameter must be large enough to allow the filter to react to changes in the 

conditions in the upper atmosphere. 4) There should be enough process noise driving the 
clock frequency errors to model the oscillator Allan variance. 
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Analytically, the process noise for the SA stales is: QD sa = block diag(JT j 



where the blocks are initialized to the steady state SA covariance matrix if the satellite is 
new in track, or the accumulated process noise appropriate for an interval of length 
A = t k+l - t k if it is not a new satellite: 



'y if v' h rowofSl(t„.„t t )= 0 
^ otherwise 



(93) 



where the steady-state covariance is: 



o 



sa_rng 



(94) 



a sa_rng is the ste ady state standard deviation of the sa process (currently 75 m), 
a saj»ei ~ a sa_m g m n ls the velocity sigma, and & n is the natural frequency of the sa 
10 markov (currently 300 seconds). The covariance for the existing satellites is: 
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c. b 
b a 

-2<re7„A 



C = 



1 e_ 



(2 a 7„„A)+-^ i sin(2^A) 
(2^A) + ^sin(2sr n ,A) 



2 I 
nd/ar n 



1 

— + COSl 



1 



- COS 



27 



(95) 



The process noise for the carrier cycle ambiguities is zero unless the satellite is new 
or a break in track has occurred, in which case it is set to a very large value. Thus: 



5 where: 



Pv = 



QD pdr = Le9 -diag(pj) 

1 if v th , v th element of Sl(/* +1 , ^ ) = 0 or /om_ tod v = 1 
0 othei~\vise 



..,(96) 



The GNS state can be initialized or reset by any one of three different ways: 

1) The GNS state can self-initialize (or reinitialize) using two GPS point solutions. 

2) Initially, the state can be set using preloaded orbital insertion data. 

10 3) The state can be set (or reset) by a command upload of a NORAD state vector. 
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Two GPS point solutions are used to initialize the state and covariance partitions for 
all of the dynamic states except for drag. In each point solution, the system solves for the 
receiver position and clock parameters. It assumes constant velocity, the two point solutions 
are therefor treated as data for estimating the position, velocity, receiver clock error, and 
5 receiver fractional frequency error at the second epoch. At least 4 satellites must be in track 
before point solutions can be obtained which are defined completely by the data. With the 
inclusion of prior information on the clock and fractional frequency error, the technique 
works with three satellites. The information on the clock is obtained from the constraint that 
the propagation delay is around 60 milliseconds (+/- 20 milliseconds) and is applied in the 
10 non-linear least squares point solutions. The information on the fractional frequency error is 
obtained from preflight calibration and is applied in the combining procedure. 

Under the assumption of constant velocity, and assuming there is no MIC control 

being applied, equation (9) may be approximated as T = - £ , then: 



HO ' 
*>>) 
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-I- At 
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0' 


0' 
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At 
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0' 


0' 


1 


0 



H 



v(4))' 

<*(',)) 

4h) 



(97) 



15 



z= [zpZj] withZy = 7^//)] is then calculated from the position-clock 

point solutions at epochs i = 1,2 (as described below). The point solution estimate errors 

are regarded as independent (although this is not true for the errors caused by selective 
availability), to enable the system to solve the Bayesian least-squares problem: 
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z=Hx+v 

■ x~A'(a,I ) 
v~W(0,R) 



where H is as above, JL, anci ^ represent any prior information on X , and R is the 

measurement noise associated with uncertainty of the single point solutions given below. 
The solution to this problem is obtained using the same method that was used to perform the 
5 measurement update, namely equations (67) and (68), This is a square root procedure 

Z^n 1/2 ^ 1/2 
= /, 2 i anc * 

R = R 1/2 R 1/2 where and R l/2 are upper triangular. Since only prior information 

and // (partitioned consistent with the partitioning of x in equation (97)) are of the form: 
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00 I 
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0' 


0' 


oo 
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0' 


0' 


0 


a 






three 


satellite solutions 


oo 


four 


or 


more satellite 



a = < 

//p[o',o\o,*]' 

three satellite solutions 

b = 

0 four or more satellite solutions 



where £ Q and a £ are the prior estimate of the fractional frequency error and its associated 
one-sigma uncertainty. The noise factor is given by: 

Y/ /2 o 



R 1/2 = 



0 



r l/2 
1 2 



5 where T^ 2 are the covariance matrices from the point solution nonlinear least-squares 

procedure described below. The aposteriori covariance and state estimate are nexed mapped 
into the Kalman filter order. 

The point solutions used above are arrived at by the following a nonlinear least 
squares solution to the point navigation equations. 
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Starting with an initial guess for the receiver position in the inertial frame (?) and 
the receiver clock error (jf) , the initial guess for the clock bias is obtained by subtracting 

the estimate of the transmit time corrected for a mean propagation time of 60 milliseconds 
with the receiver reference time from the tracking processor: 

f =t k -(tl + 0.060) 

The true time at the receiver r = t+ T isthen determined: 
For tracking channels v = 1, /7 max , the implicit light time equation is solved: 



A v = iLv( r v). f(r) | 



c 1 

t v =t-A v (98) 
r v (r) = T;;(r)r c ;(r) 



10 where r^ s (r) is the satellite position according to the message ephemeris 

After convergence, the unit vectors are computed: 

_ r'(r-A')-Kr) 
" -H-A 2 )-f(r)| ( 

The sensitivity matrix is then created: 
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M = 



-(u ""»)' 



(100) 



If channel v is not tracking a satellite, that row of M is zeroed out. 
The measurement vector is formed: 



y = [^ 2 ,-"X, m „>o] 

where: y v = p - c[a v + (f - r(r v ))] 



(101) 



predicted pseudorangc 



where /? * is the raw pseudo-range measurement and T{t v ) is the satellite clock bias 

computed from the message data. (The SA and the ionosphere are ignored in the foregoing 
initialization procedure.) 

Next, errors in the state are determined and the state is updated: 



10 



f 




f 


T 

- - 


new 


f 



+ By 



(102) 



old 



where B is the Kalman gain associated with the Bayesian least-squares problem: 
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y = Mx + v 
xTN(0,D) 
v~A>(0,E) 
ool 0 
0' e 



a r three satellite solution 

od four or more satellite solution 



where E^ 2 = diag(VKRJE) is a diagonal matrix representing the average user effective 
ranging error and <j t = .0.020 represents the uncertainty in the coarse estimate of clock 

error obtained from the constraint that the range delay is aproximatly 60 milliseconds. (This 
5 is dominated by SA if it is on.) The above process is repeated until the estimates are not 

changing rapidly, i.e. |fiy| < test . At convergence, the error covariance factor using the 

algorithm of equation (67) is selected as P 1/2 . This is used to form the measurement noise 
factor for the combining step described previously. 

The Delta guidance system will place TIMED at a position and velocity specified in 
10 the Conventional Terrestrial System (CTS). 
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C'"(>) 



V 



(/•) 



represents the true position and earth-relative velocity of the 



CIS 



6x1 



TIMED CM in the CTS system and W^"(/) represents its nominal value. Then the orbital 
insertion error appearing in the Delta specification will be of the form: 



5 where r ins is the true orbital insertion time, defined to be the time immediately after all 
thrusting (including pyrotechnic separation devices) fs over, and t™ m is its nominal value. 

By specifying the insertion error in this way, there is no penalty for placing TIMED in the 
correct position and velocity at a different time than intended. The desired orbital insertion 

vector w™ m (/™") is a point fixed in the Conventional Terrestrial System (CTS). If the 
10 launch time slips, this aim-point will remain constant, implying that the right ascension of 
the ascending node of the TIMED orbit will rotate at Earth rate. The six-vector w nom ( t nom ) 

els \ ins J 

will be loaded into memory on the spacecraft prior to launch. The actual orbital insertion 

time ( T in S )> which may be substantially different than the nominal time (/™ m ), is 

communicated to the TIMED GNS via a measurement of the TIMED disconnect switch time 

15 ( r dsu )• The difference = t dsw - r ins is the error in time tagging the insertion event and 

is assumed to be small (at most a few seconds). The TIMED Guidance and Navigation 
System (GNS) position and velocity state vector is defined in the Conventional Inertial 

System (CIS). Therefore w™ m (/™") must be transformed into the CIS using a 
transformation valid at / . . : 




(103) 



20 
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where: 



A(/) = 



TS(/) 
0 



0 



I 



o 



The transformation from Earth-relative to inertial velocity is built-in to A(t) . The 
synchronous error of the TIMED position and velocity state at /^ M , is then: 

= w c/J (^)+ ^(rJ-(r-(/^) - f te )- A(^)-wr(c) 
- [A(^) + A(^„) • (^,. - r J] • wr(C') 

= w c ,(r ;m )- A(r,)w;:(c)+ w^rj.^ - t^)) 

-A(r (/M )-wr(cr)-7;t' 



5 where: 



A(r) = 



cwk x l 


0 


I 0" 


0 




>,*} I. 



bias dsw ins 
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^ w c/-j 5 nas lo °e written in terms of the orbital insertion error defined in equation 
(103), since w eif (r te )= A(r in J- w c/J (r.„J : 

+ - r(t dsw j)- A(^s)<T(C m )- r b t 

= A(rJ.*w£ + w rt (rJ-(2£T - jnfcj) 

-A(rJ'w;:(c)'f 
= A(r iu ).K;-w cir (rJ.r(g 



where W c/J (/ Aw ) = f (w c/j )(f^ u> ) is the model predicted rate of change of inertial position 

5 and velocity (largely driven by local gravity). Arranged in matrix form to clearly show the 
error source contributors as: 



2tO 



A(rJ 
0 

~ A ('-J 
0 



w(r, I )-A(r,,).w:r(c)' 



0 



iJ-4).w;;(c) 



bias 

5\\ 



rrtdSM' 



0 



A/ 



cts 



nndsw 
1 bias 
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The following approximations have been used 

A(0= A.(r lnf )+ A(r ins )(t dsw - r ms ) 
= A(r /Bf )+ A(r /W )7^T 



to substitute the computable A/ x , w, v and A/ x for the corresponding quantities 

evaluated at (the unknown) r frtf , realizing that terms like 7^A(r. /w )^w^ are higher 

5 order. The covariance structure for the synchronous errors consistent with the initialization 
in equation (104) is then: 

(105) 



cov(^wf-) cov(^wf~, zj 
[symmetric) cov{T) 



= M 



covfcXZ) 0 0 

0 cov(r) cov(r,r b t:) M 

0 (symmetric) cov(2^*) 



The actual correlation between the errors in the disconnect switch time tagging of the 
10 insertion time and the receiver clock error (cov(t\ 7^')) depends on whether and how 

recently the C&DH clock has been synchronized with the GNS lpps. It is assumed to be 
zero. 

In addition to the initial setting of the position and velocity state vector, the system can reset 
the state vector by command based on NORAD radar tracking. This command consist of 

15 CIS position and velocity estimates along with a defined time of validity {j up}< j) - 
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The firsi sicp in implementing this upload iu Lraiiaibrm ilic orbital cicnicnls 
provided by NORAD on its two-line element set format into the CIS position and velocity: 



, _ ffaj (106) 

where — 

d6 



where 6 are the orbit-defining quantities in the two-line element set, and \ T upId ^j is the 

5 nonlinear transformation implemented in SPG4 that takes these quantities into the inertia] 
position and velocity at time v ld . The derivative of the transformation E can be found 

analytically or by the numerical differentiation of SGP4. The error Sw n °™ d {j up}d ^ is not 

identified as synchronous or asynchronous because there is no error associated with r ld 

since it is a defined time at which the SGP4 transformation is evaluated. The error (if any) 
10 associated with time-tagging the NORAD epoch will propagate through the transformation 

^(•) . In this case, the synchronous and asynchronous errors are equal, so a generic 
designator is used. 

Since the Kalman filter time update propagates between receiver times, the time of 
validity must be transformed from a true time to a GNS receiver time, using the process 
15 illustrated below: 
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where r "'(•) is the inverse of the function that takes GNS reference times into their 
corresponding true times. Specifically, the GNS time corresponding to T upld is calculated 
from the relation: 

/ _ t gns 

l k l upld 

= 1 + s k 

T k ~ T upld 

where: 

t k = current receiver time (know) 

r k = true time corresponding to t k 
T u P id = timQ °f validity of uploaded state acording to NORAD 
*upM = receiver clock reading coiTesponding to r upld 

£ k = current estimate of receiver frequency bias 



which is solved for t*^ d to obtain: 



= t k -{T k -f k -T upld ){\+S k ) 



(107) 



10 The foregoing completely defines the TIMED state upload. The position and 

velocity are given by: 
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w Z ns (t8" s ) = W norad (r ) 
W cis \ 1 upld ) ™ cis \ 1 upld ) 

and the start time for the time update (t las( ) is set to t*^ d from equation (107). 

The covariance matrix implied by this initialization process is derived s follows. The 
5 synchronous error of the GNS state is: 

~ W cis ( T upld ) ~ W cis ( T upld ) 

= ^r(^) 



10 



indicating that the GNS inherits the synchronous error of the NORAD uplinked state vector, 
which is the generic error defined in equation (106). The appropriate transformation of the 
error state is then: 



(t gns V 

\ L upld j 



I 0 

0 1 



6w nc 

as 



norad 



The appropriate transformation of covariance is: 



it 



cov^w) 
symmetric) 



cov{S s w,T) 
cov(T) 



l upld 



cov(*C-) 
(symmetric) 



0 

cov(r) 
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The foregoing defines only the GNS position, velocity, and clock error covahancc. 
Default values are used for the rest of the state space. 

Short Term Propagator 

5 The short term propagator creates the navigation data that is placed on the TIMED 

data bus every second. This short term data includes (in the order of calculation): 

• Data Validity Flags (5), 

• CIS frame Position and Velocity, 

• Time, 

10 • Sun Vector, 

• Geodetic Latitude, Longitude, Altitude, 

• East, North, Up Velocity, 

• Day/Night Event Notification Flag, 

• South Atlantic Anomaly (SAA) Event Notification Flag, 
15 • Polar Region Event Notification Flag, 

• Ground Station Event Notification Flags (GSENF) 

The methodology used to create these data is described by the following derivations wherein 
reference is made to the dynamic propagation rutine, a primary part of the Kalman Filter 
(KF) that performs the time update for the KK Crank. This routine is capable of operating in 

20 both covariance propagation and non-covariance propagation modes. 

During the short term propagation, the dynamic propagation routine operates in the 
noncovariance propagation mode. 

These data are computed after the time update in each KF Crank interval. This is part 
of the "Cleanup" of the KF Crank interval. During each KF Crank, the short term data is 

25 computed for a 180 second interval beginning at the start time of the next KF Crank. Thus, 
the basic flow of the KF Crank for a single crank interval is: 

• KF Crank Start at T f 

• Command Input 

• Measurement Update, 
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• Time Update (TU) (slate propagated to 7^ lR0 ) 

• Short Tcmi Propagate (for 7; +lg0 to T i+360 ) 

In the following description, the short teim data is valid for the period T i+J80 to T i+360 . 
The short term data cannot be created for the same interval as the KF Crank because the KF 
5 Crank takes longer than 1 second to complete. Therefore, the short term data for the for the 
first several seconds of the current interval is not available when it is needed. 

The short term propagator begins with the propagated dynamic state vector X" 
(position(3), velocity(3), drag, clock, frequency), state vector derivative X~ , and covariance 
factor matrix LT . These variables are created by the time update part of the KF Crank and 

10 are valid at time T i+/80 . To simplify the notation, the superscript denoting "before 
measurement" is dropped. Assuming T s0 = T i+/80 the required calculations are: 

• Compute Data Validity Flags 

• Propagate State to T^ I80 and T l+360 

• Compute Polynomial Fit Coefficients 
15 • First 180 Second Loop: at every second 

• Compute X(r) 



• Compute 




20 



• Compute Sun Vector 



• Compute OTS Frame Position, Velocity, and Sun Vector 



25 



• Take Q x R out of velocity in X 



CIS 



110 



3NSDOCID: <WO 00636*6A1J_> 



WO 00/63646 



PCT/US00/10197 



• Compute Geodetic Latitude, Longitude, and Altitude 

• Compute East, North, Up Velocity 

5 • Compute Day/Night Event Notification Flag 

• Compute South Atlantic Anomaly Event Notification Flag 

• Compute Polar Region Event Notification Flag 

10 

• Compute Ground Station Event Notification Flag (GSENF) 

• Second 180 Second Loop: at every second 

• Compute X(/) 

15 • Compute T™(t) 

• Compute GSENF 

• Shift GSENFs 

20 The data validity flags are computed using the L matrix. They are constant 

throughout the 180 second interval. The following values are used to compute the one 
sigma values: 

The covariance matrix P = L L r and that L is lower triangular. This means that 
the diagonal values of P are given by P(i 9 i) = ^ - This y ields the one sigma 
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values for position, velocity and time. The data validity flags are then set from these one 
sigma values according to the following: 

• if all three values are < the position validity value, the position and event 
validity flags are set to VALID, 

5 

• if all three 0" velocity values are < the velocity validity value, the velocity valid velocity flag 
is set to VALID, 

• if the a chck value is < the time validity value, the time and sun vector validity flags are 
10 set to VALID, 

The same dynamic propagator that is used in the Time Update of the KP Crank is used to 
propagate the state to T i+I80 and T i+360 .. (See equation (34) and the derivations leading 
thereto.) This gives state vectors and state derivative vectors at SO, S0+180(S1), and 
15 S0+36.0(S2 ). At. this point the system has the following information: 

X(S0)' 9 X(S0), X(Sl), X(S1), X(S2) 9 X(S2) 

The system requires short term data at a 1 -Hertz rate, therefor the data is interpolated 
between SO, SI , and S2. This is accomplished by fitting a 5th order polynomial to the curve 
described by: 

20 X(S0) 9 X(S0) 9 X(Sl) 9 X(Sl) 9 X(S2), X(S2) . 

This structure is used to interpolate data for the state vector and the state vector 
derivative defined as follows: 
For the state vector: 

X(/) = A+B-f + C-r +D-/ 3 +E-f 4 + F-/ 5 (108) 
25 For the state vector derivative: 

X(/) =B+2-C-f +3-D-/ 2 +4-E-/ 3 +5-F-/ 4 (109) 
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The interval between SO and SI is the same as the interval between SI and S2 and 
denoted by /. By letting the time(t) for the above equations go from -180 to +180 instead of 
0 to 360, the system can then use the following six equations to determine the six 
coefficients: 

X(t = -/)=A-B/+C/ 2 -D/ 3 + E/ 4 -F/ 5 
X(t = -/) = B- 2C/ + 3D-/ 2 -4E/ 3 + 5F/ 4 

X(/ = 0) = A 
X(t = 0)=B 

X(/= /) = a + B-/ + C-/ 2 + D-/ 3 + E-/ 4 + F-/ 5 
X(t = I) = B+ 2-C-J + 3-D- 1 2 + 4-E-/ 3 + 5 F / 4 



This yields the following equations for the six coefficients: 
A = X(S1) 
B = X(S\) 

C = (l// 2 ) • (X(S2) - 2 ■ X(Sl) + X(S0) - 1 ■ X(S2)/4 + / • X(S0)/4) 
D = (l/(4 • 7 3 )) • (5 - X(S2) - 5 • X(S0) - 8 • / • X(S\) - I ■ X(S2) ±J, X(S0)) 
E= -(l/2-/ 4 )-(x( 1 S2)- 2-X(Sl)+ X(S0) - I ■ X(S2)/2 + I ■ X(S0)/2) 
F = (l .0/(4 • J 5 )) ■ (3 • X(S2) - 3 • X(S0) - 4 - / • XGSl) - / ■ X(S2) + I ■ X(S0)) 
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The foregoing provides the coefficients required to compute X(/) at any point by 
looping through the first 180 second interval and proceeding as follows at each of the 1 80 
(*so+i to *so+i8o) Points. 

Using the previously described points, equation (108), 

5 X(t) = A+B-f + C-r 2 +D-/ 3 + E-/ 4 + F-/ 5 , X(/) is computed each time ( t) 
goes from -180 to 180 while the trajectory time goes from SO to S0+1 80. 

The same routine is used to compute that is used in the Dynamic Propagation 
Routine (see also equations (71) and (72) and the comments therefor). 

The sun vector, S CIJ , in conventional inertial System, CIS, coordinates is computed 
10 using the same routine that is used in the Dynamic Propagation Routine. 

With and S c/J computed, the system determines position, velocity and sun 

vector in the conventional terrestrial system, CTS, frame: 

The position vector in the CTS frame is determined by: 

P«, =T c c ;;X cis {position) (110) 

1 5 The velocity vector in the CTS frame is determined by transforming the coordinate 

system from CIS to CTS and removing the earths relative velocity from the inertial velocity. 
The velocity vector in the CTS frame is: 

V« = • Xjvelocity) - Ci earlh<cts x P c „ (ill) 



The matrix multiplication is eliminated because: 







0 


- ® earth 


o" 






20 


CI V P - 

earth, as as 


®eanh 


0 


0 




Pi 






0 


0 


0 




A. 



which simplifies the process and results in: 
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(n^ tt x Pj(0 = -^-P„(2) 

I ^ earth cis X Pc«)(2)= ®ear,„ KM 



The sun vector in the CTS frame is: 



C/5 - l dS '^CIS 



(112) 



The Geodetic Latitude {lat) and Altitude are computed using a routine that is the 
5 same as used in the Dynamic Propagation Routine. The Longitude (Jon) is computed using 
the standard C library routine atan2: 

longitude = a tan 2(P C/J (y_ axis), P cfs (x_ oxw) (113) 

To compute the east, north, and up velocity, the system computes the transformation 
from CTS coordinates to East, Down, and North @ Equator (EDN) coordinates: 



10 



nnedn __ 
*cts ~~ 



cos{lon + */Q sm{lon + 0 

- sin {ion + co^lon + *fy 0 

0 0 1 



It then computes the transformation from EDN to East, North, Up (ENU) coordinates: 



npenu 



0 0 
0 cos(^-/o/) sin(^-/*rf) 

0 - sin(^ - lat) cos{Y 2 - lat) 



npenu np 
* edn * cis 



V = T enu • V 

enu edn cis 
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The system uses the Geodetic Latitude, therefor "UP" is in the local vertical direction. 

Because the normal to the surface in the ENU coordinate frame is given by the vector 
[0,0,1], the system transforms the vector to the CTS frame using: 



U , = T cts 

cts enu 



which provides: 



U„(l)-T~(l,3)=TT(3,l) 

U M (2)=T«(2,3)=Tr(3.2) (115) 

u ctt (i)=T;;„(3,3)=T;r(3,3) 



This provides a vector normal to the surface of the earth in the CTS coordinate frame. The 
system then creates: D = U cts • S cls . If U cts • S cls > 0.0, then the latitude and 

longitude point is in sunlight and the flag = 1 . When U cls • S cu < 0.0 , the flag = 0. 

The South Atlantic Anomaly, SAA is described in the GNS by 12 longitude and 
10 latitude points which are vertices, in clockwise order (looking down at the earth), of a 
convex polygon inside which lies the SAA. When these values are loaded into the KF 
database, the system computes vectors for each side using S- = V. + ] - V. , where 

V. = {longitude., latitude^) . These are longitude and latitude points (NOT latitude, an 

longitude) and are ready for use in the Short Term Propagation loop. 
15 To determine if the spacecraft is within the SAA polygon, the system performs the 

following operations for each side of the SAA after setting the SAA event notification flag 
to 1 EVENT'. 
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The spacecraft vector, S sc , , is computed using the previously computed spacecraft 

latitude and longitude {Loh sc , Lat sc ) , which is in the form of the vector from V. to the 
spacecraft location: 

Ssc^iLon.^Lat^-V, (116) 

A z-coordinate (=0) is appended to each vector and used to successively compute the 
cross-product of a side of the SAA polygon with the spacecraft vector: 





0 




0 


0 








Q = s, x s 5C ,. = 


0 




0 


0 






(117) 




a. 




r s >? 




0 


0 





The system uses the right hand rule, i.e., if Q z > 0 the spacecraft lies 
outside the SAA and the SAA Event Notification flag is set to 'NO-EVENT'. The system 
10 then exits the loop. 

The Compute Polar Region Event Notification flag is set using a simple comparison, 

i.e., if \Latitude\ > Value , the flag is set to 1 . 

The location of each ground station is described by a latitude and longitude. When 
these values are loaded into the KF Database, the system computes CTS coordinates 
15 V C5| . as for each of the ground stations using the ground station latitude and longitude and 

an altitude 0. It then computes the CTS vector N CSi ; cts normal to the surface of the earth at 

the ground station coordinates using the following steps: 

Multiplying the two matrices of equation (114) reduces the operation to: 
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ft 

0 = —-lat 

2 

71 

1 = Ion + — 
2 

cos(^) sin(^) 0 



*cts 



- cos(#) sin(^) cos(#) cos(^) sin(#) 
sin(#) sin(^) - sin(#) cos(^) cos(#) 



The normal to the surface in the ENU coordinate frame is given by the vector [0,0,1] 
as established by equation (115). The vector is transformed to the CTS frame 



N = T cts 



This yields: 



N G *, C/ ,(1)= T-(U)= TT(3,l)= sinO?)sinM 

N C5l , ctt (2) = T-(2,3) = Tr(3,2) = - sin(tf)cosM (l IS) 

N G5 , ctt (3) = T-(3,3) = 1T(3,3) = cos(0) 



The values are now ready for use in the Short Term Propagation loop. To determine 
if the spacecraft is above the horizon relative to each ground station, the system computes 
the dot product of the vector from the ground station to the spacecraft and the ground station 

normal vector: D= (P , - V^ c . )• A^ c 

\ as GSi,cts J J GSi.cts 

10 If D > 0 , the spacecraft is at or above the horizon and the GSENF is 1 . 
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The system computes the second 180 second loop using the previously described 
operations, i.e.: 

X(/) is computed each time ( / ) goes from -180 to 180 while the trajectory' time 
goes from SO to S0+180. 
5 X(f)=A+B-/+ C-t 2 + D/ 3 +E-/ 4 

The same routine is used to compute that is used in the first 180 second loop. 

The same routine, equation (118), that is used in the first 180 second loop is used to 
compute the GSENF. 

The system then shifts the GSENF as the short term propagation continues to provide 
10 GSENFs for each ground station for a period of 360 seconds. The extra 180 seconds allows 
the system to predict rise or set times up to 1 80 seconds in advance of the time it 
outputs the data. The GSENFs which are output may need to be set before the spacecraft 
rises above the horizon and may need to be unset before the spacecraft has set below the 
horizon to satisfy Rise/Set Lead Time requirements. The system predicts rise time so that 
15 instruments on board the satellite can be ready to send data when the satellite rises above the 
horizon instead of spending time gaining acquisition to the ground station. The system also 
computes set lead time so that the satellite may know when to start shutting down the data 
transmission sequence. The time shifts for the Rise Lead Time R h and Set Lead Time S h 

are GNS parameters. 
20 For example, if R lt = 4 seconds and the GSENFs are 

00000011111. ..Ill 
the time shifted GSENFs are 

0 0 01 1 1 1 1 1 1 1 ... 1 1 1 

Or if R lt - 5 and the GSENFs are 

25 1 1 1 1 ... 1 1 1 1 1 1 1 1 0 0 

the time shifted GSENFs are 

1111. ..1110000000 
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There are several situations that can occur in the stream of GSENFs during the 360 
second interval. 

1 ) There are no GSENFs set. These flags do not need to be modified for Rise/Set Lead 
Times. 

5 2) The GSENF stream started unset, changed to set and stayed set for the remainder of the 
interval. These flags need to be modified as in the first example above. 

3) The GSENF stream started set, changed to unset aud stayed unset for the remainder of 
the interval. These flags need to be modified as in the second example above. 

4) The GSENF stream started unset, changed to set, and changed back to unset. These flags 
10 need to be modified, first, as in the first example, and next, as in the second example 2. The 

system performs the modifications in that order to cover the situation in which the GSENFs 
are set for a period of time < S it - R n . 

The foregoing is considered as illustrative only of the principles of the invention. 
Since numerous modifications and changes will readily occur to those skilled in the art, it is 
15 not desired to limit the invention to the exact construction and applications shown and 
described. Accordingly, all suitable modifications and equivalents may be resorted to and 
are considered to fall within the scope of the invention and the appended claims and their 
-equivalents. 

What is claimed is: 
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Claim 1. A satellite navigation system, comprising: 

a radio frequency receiver incorporating a plurality of channels for 
concurrently processing code containing signals received from a plurality of 
global positioning system satellites; and 

a microprocessor for executing software which manipulates data processed by 
said plurality of channels according to an extended Kalman filter to 
autonomously generate estimates of position, velocity and time. 

Claim 2. A satellite navigation system as defined by claim 1, wherein said radio 
frequency receiver comprises: 
a plurality of antennas; 

a plurality of preamplifiers, driven by each one of said plurality of antennas; 
a separate downconverter for each one of said preamplifiers for processing 
the amplified signals produced thereby; and 

an application specific integrated circuit for producing tracking data for each 
of said global positioning system satellites from said amplified signals 
processed by said downconverters. 

Claim 3. A satellite navigation system as defined by claim 2, wherein each of said 
downconverters comprise: 
a phase locked loop frequency synthesizer; 

a reference oscillator for driving said phase locked loop frequency 
synthesizer; 

a mixing means for combining an output of said phase locked loop frequency 
synthesizer with an output of said preamplifier; 

an intermediate frequency band pass filter for processing the output of said 
mixing means; and 
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a quantizer for producing a pair of complementary signals from an output of 
said intermediate frequency band pass filter. 

Claim 4. A satellite navigation system as defined by claim 3, wherein said application 
specific integrated circuit comprises: 

a latch circuit for sampling said pair of complementary signals from two of 
said downconverters and producing therefrom a pair of complementary 
quadrature signal samples and a pair of complementary non-quadrature signal 
samples; 

a plurality of tracker channels; and 

a data router for distributing each of said complimentary signal sample pairss 
to a different one of said tracker channels. 

Claim 5. A satellite navigation system as defined by claim 4, wherein said tracker 
channels each comprises: 
a plurality of accumulators; 

a separate latch register for outputing the contents of each of said 
accumulators for outputting the contents thereof; 
a coarse acquisition code generator; and 

a differential mixer for combining said signal samples with a code created by 
said course acquisition code generator to lock the phase and frequency of said 
code carried by said signal samples with said code created by said course 
acquisition code generator. 

Claim 6. A satellite navigation system as defined by claim 1, comprising: 

an application-specific integrated circuit for downconverting code containing 
GPS signals received from said plurality of global positioning system 
satellites and tracking said code containing GPS signals. 
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Claim 7. A satellite navigation system as defined by claim 6, wherein said 
microprocessor comprises: 

a tracking processor for controlling said application-specific integrated circuit 
GPS code tracking functions; and 

a navigation processor for implementing said software execution of said 
extended Kalman filter. 

Claim 8. A satellite navigation system as defined by claim 6, comprising: 
a plurality of input channels for producing GPS signal samples; 
a plurality of tracker channels for accumulating code carried by said GPS 
signal samples ; 

a data router for distributing said GPS signal samples to said tracker channels 
as a function of tracker channel availability and said global positioning 
system satellite originating said GPS signal samples; 
a timing means synchronizing the circuits of said application-specific 
integrated circuit; and 

a bus interface for coupling said accumulated code to said microprocessor 
and for receiving controlling data from said microprocessor. 

Claim 9. A satellite navigation system as defined by claim 8, wherein each of said 
input channels comprises: 

a first exclusive NOR gate for processing a pair of complementary 
non-quadrature signals derived from one of said code containing signals; 
a second exclusive NOR gate for processing a pair of complementary 
quadrature signals derived from one of said code containing signals; 
a first multiplexer for producing a non-quadrature signal sample in response 
to an input from said first exclusive NOR gate; 

a second multiplexer for producing a first quadrature signal sample in 
response to an input from said second exclusive NOR gate; 
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a third multiplexer for producing a second quadrature signal sample 
complementary to said first quadrature signal sample from one of said 
complementary non-quadrature signals and one of said complementary 
quadrature signals; and 

an output means for a non-quadrature signal sample complementary to said 
signal sample produced by said first multiplexer from said one of said 
complementary non-quadrature signals. 

Claim 10. A satellite navigation system as defined by claim 9, wherein each of said 
tracker channels comprises: 
a numerically control oscillator; 

means for combining a cosine function of said numerically control oscillator 
with said non-quadrature signal samples received from said IF input channel; 
means for combining a cosine function of said numerically control oscillator 
with said quadrature signal samples received from said IF input channel; 
means for combining a sine function of said numerically control oscillator 
with said non-quadrature signal samples received from said IF input channel; 
means for combining a sine function of said numerically control oscillator 
with said quadrature signal samples received from said IF input channel; 
means for summing a product of said cosine function combined 
non-quadrature signal samples and said sine function combined quadrature 
signal samples; 

means for summing a product of said cosine function combined quadrature 
signal samples and said sine function combined non-quadrature signal 
samples; 

a non-quadrature flip-flop for producing a plurality of parallel output signals 
of different values from an output of said means for summing a product of 
said cosine function combined non-quadrature signal samples and said sine 
function combined quadrature signal samples; 
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a quadrature flip-flop for producing a plurality of parallel outputs signals of 
different values from an output of said means for summing a product of said 
cosine function combined quadrature signal samples and said sine function 
combined non-quadrature signal samples; 

a coarse acquisition coated generator; a coarse acquisition code numerically 
control oscillator for controlling said coarse acquisition code generator; 
a multibit shift register responsive to outputs from said coarse acquisition 
code numerically control oscillator and said coarse acquisition code generator 
for producing a plurality of time shifted coarse acquisition codes for 
combining with said plurality of outputs from said non-quadrature and 
quadrature differentiators; 

a plurality of accumulators, each of which accumulates code resulted from 
said combination of said time shifted coarse acquisition codes and said 
plurality of outputs from said non-quadrature and quadrature flip-flop; and 
a latch register for each of said accumulators for setting and output thereof. 

Claim 11. A method for orbital platform navigation, including the steps of: 

sensing the data transmitted by a plurality of global positioning system 
satellites; 

determining from said sensed transmitted data through the application of a 
software implemented Kalman filter, data defining current time and position 
and velocity of said orbital platform; and 

initiating and controlling the operation of systems on-board said orbital 
platform based on said data defining current time and position and velocity. 

Claim 12. A method for orbital platform navigation as defined by claim 1 1, wherein 
said Kalman filter utilizes previously determined data defining time and 
position and velocity when performing said step of determining from said 
sensed transmitted data, current time and position and velocity. 
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Claim 13. A method for orbital platform navigation as defined by claim 12, wherein 

said step of determining current time and position and velocity includes steps 
executed by said Kalman filter of calculating a square-root covariance form. 

Claim 14. A method for orbital platform navigation as defined by claim 13, wherein 

said Kalman filter preforms the steps of determining the measurement of slant 
range from said orbital platform to each of said global positioning system 
satellites providing data used in determining said data defining current time 
and position and velocity. 



Claim 15. A method for orbital platform navigation as defined by claim 14, wherein 
said plurality of global positioning system satellites comprise at least three 
satellites. 

Claim 16. A method for determining orbital platform navigation as defined by claim 15, 
wherein said Kalman filter preforms the steps of providing a plurality of state 
space determinations arranged in first and second partitions; 
said first partition including states which directly determined the orbital 
platform position, velocity, drag parameter error, receiver clock and 
frequency bias; and 

said second partition includes states not directly affecting position and time. 



Claim 17. A method for orbital platform navigation as defined by claim 1 6, wherein 
said steps of providing a plurality of state space determinations include the 
steps of: 

calculating a position state space; 
calculating a velocity state space; 
calculating a drag parameter error state space; 
calculating a clock error state space; 
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calculating a frequency error state space; 

calculating a selective availability and Markov process state space; and 
calculating an integrated carrier phase bias state space. 

Claim 18. A method for orbital platform navigation as defined by claim 16, wherein 
said steps of providing a plurality of state space determinations include the 
steps of: 

calculating three position state spaces; 
calculating three velocity state spaces; 
calculating a drag parameter error state space; 
calculating a clock error state space; 
calculating a frequency error state space; 

calculating 24 selected availability and Markov process state spaces; and 
calculating 12 integrated carrier phase bias state spaces. 

Claim 19. A method for orbital platform navigation as defined by claim 16, wherein 
said steps of providing a plurality of state space determinations include the 
steps of: 

Calculating a plurality of dynamic states including, a position state space, a 
velocity state space, a drag parameter error state space, a clock error state 
space, and a frequency error state space; and 
dynamically coupling said dynamic states to each other. 



Claim 20. A method for orbital platform navigation as defined by claim 13, wherein 
said Kalman filter preforms said step of updating said data defining current 
time by executing steps propagating the state and covariance between 
measurements. 

Claim 21. A method for orbital platform navigation as defined by claim 1 9, wherein 
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said step of propagating the state uses lime as an independent variable. 

Claim 22. A method for orbital platform navigation as defined by claim 13, wherein 
said steps of determining current time and position and velocity include 
update steps executed by said Kalman filter comprised of propagating state 
and covariance from the present to the next anticipated measurement for the 
dynamic states of position, velocity, drag, and clock, each of which is 

expressed as a function x dyn = f(x dyn ,/) where t equals time, according to 
the nonlinear function F d — where the covariance is F . 

" ^dyn 



Claim 23. A method for orbital platform navigation as defined by claim 22, wherein the 
deterministic part of said covariance factor F used in preforming said update 

steps is defined by P£ ]/Jc = . , t k )' , where P^, is the state 

error covariance after a deterministic update and prior to a stochastic update. 

Claim 24. A method for orbital platform navigation as defined by claim 23, including 

the Kalman filter executed steps of propagating state and covariance from the 
present to the next anticipated measurement for the states and biases which 
obey linear dynamics according to the function 
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Claim 25. A method for orbital platform navigation as defined by claim 24, wherein 

said states and biases which obey linear dynamics include availability states 
for said global positioning system satellites and phase biases. 
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